program main 



main program.txt 



integer, parameter :: maxres=1000 
integer, parameter :: mxratm=200 
integer, parameter :: maxsize=200 

character (len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

Theory='HF/3-'2lG' 

level =1 

cut=2 

cal 1 MFCC (Theory , 1 evel , cut , maxres , mxratm , maxsi ze) 
end 

c=============================================================== 

subrouti ne mfcc (Theory , 1 evel , cut , maxres , mxratm , maxsi ze) 

c ^ — Molecular Fractionation with conjugate caps c 



c c 

c coord: cartesian coordinates of each atom in each residue c 

c res_atomnum: the number of atoms in each residue c 

c res_atomname : the name of the atom in each residue c 

c res_name: residue name c 

c charge: charge of each residue (only 5 amino acids can have c 

c a charge and we just double-check them c 

c endcharge: check the charge of the terminals (like, nh3 at c 

c N-terminal means positive charge and COO at c 

c C-terminal means negative charge c 

c numres: the total number of residues in the protein c 

c Theory: method in the ab initio calculation (like, HF, b3lyp...)c 

c Basisset: basis set used in the ab initio calculation c 

c cut=0 cut between CA and N c 

c levell: nh2 CH3 c 

c level 2: NH2 CH2R c 

c level 3: NHCOH ch2r c 

c level 4: nh2 CHRC0NH2 c 

c levels: nhcoh CHRC0NH2 c 

c cut=l cut between CA and C c 

c levell: CH3 --- C0NH2 c 

c level 2: CH2R -— C0NH2 c 

c level 3: CHRNH2 — - C0NH2 c 

c level 4: chrnhcoh --- conh2 c 

c cut=2 cut CONH (planar Frag) c 

c levell: CH3NH C0CH3 c 

c level 2: CH2RNH --- C0CH3 c 

c level 3: CH3NH --- C0CH2R c 

C level 4: CH2RNH C0CH2R C 

c error: if error=l, we do three-body correction c 

c if error=0, no correction c 



^* ***************** Tfr:!f************************* 
^-***********************************************************************^ 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 
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character(len=4) re5_atoniname(maxres,mxratm) 
characterO en=4) res_name(maxres) 
integer res_atomnum(maxres) 

character (len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

integer maxsize 
integer charge (maxres) 
integer endcharge(2) 

call readpdb (coord, res_atomname, res_atomnum, res_name, 
$ charge , endcharge , numres , maxres , mxratm) 

call cal energy (coord, res_atomname, res_atomnum, res^name, 
$ charge , endcharge , numres , maxres , mxratm , 

$ maxsize, Theory, level ,cut) 

call SystemCrm -f fort.* out* 100* 200*') 

end 

c 

subroutine cal energy (coord, res^atomname, res_atomnum, 
$ res_name , charge , endcharge , numres , maxres , 

$ mxratm, maxsize, Theory, level ,cut) 

^************** ******* ********************** 

^************ ************************ ***********************************^ 

c atom_symbol : the name of atoms in the small molecule c 

c ligand: the cartesian coordinates of atoms in molecule c 

c ligand_atomnum: atom number of ligand molecule c 

c Frag: the cartesian coordinates of atoms in each fragment c 

c Frag_atomname: the name of atoms in each fragment c 

c Frag^atomnum: the number of atoms in each fragment c 

c Frag_id: the polarity of each residue c 

c nopolar(l) ,polar(-l) , polar with charge(O) c 

c Cap: the cartesian coordinates of atoms in each cap c 

c Cap_atomname: the name of atoms in each cap c 

c cap_atomnum: the number of atoms in each cap c 

c Cap_id: the polarity of each cap c 

c c 

c In correction, we put two nearest residues together c 

c and treat them as one fragment (like, 12,23,34,45. . .) c 

c corr: the cartesian coordinates of atoms for each fragment c 

c corr_atomname: the name of atoms in each corr fragment c 

c Corr_atomnum: the number of atoms in corr fragment c 

c pcharge: the charge for each fragment c 

c ccharge: the charge for each cap c 

c Corr^charqe: the charge for each corr fragment c 

c p_frag: the sum of fragment energies c 

c p_cap: the sum of cap energies c 

c p_su I : the sum of sulcap energies c 

^***********************************************************************^ 

^*********** *********************************************************** *^ 

integer maxres 
integer mxratm 
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real*8 coord(maxres,nixratm, 3) 
character(1en=4) res_atomname(maxres,mxratm) 
characterOen=4) res_name(maxres) 
integer res_atomnum(maxres) 



integer maxsize 

character(len=4) atom_symbol (maxsize) 
real*8 ligand (maxsize, 3) 
integer ligand_atomnum 
integer ligand_charge 

real*8 Fraq(maxres,mxratm,3) 
character(Ten=4) Frag_atomname(maxres,mxratm) 
integer Frag_atomnumCmaxres) 
integer Frag_id(maxres) 

real*8 Cap(maxres,mxratm,3) 
character(len=4) Cap_atomname(maxres,mxratm) 
integer Cap_atomnum(maxres) 
integer Cap_id(maxres) 

rea1*8 Corr(maxres,mxratm,3) 

cha ract e r (1 en=4) Co r r_at omname (max res , mx r atm) 

integer Corr_atomnum(maxres) 

integer Corr_charge(maxres) 

real*8 Sulcap(maxres,mxratm,3) 

character (len=4) Sulcap_atomname(maxres,mxratm) 

integer SulCap_atomnum(maxres) 

integer sn 

character(len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

real*8 fgra_lig(maxsize,3) 
real*8 cgra_lig(maxsize, 3) 
real*8 sgra_lig(maxsize,3) 
real*8 gra_lig (maxsize, 3) 

integer charge(maxre5) 
integer endcharge(2) 
integer pcharge(maxres) 
integer ccharge(maxres) 
integer scharge(maxres) 

real*8 f_asymp(maxres) 

real*8 c_asymp(maxres) 

real*8 s_asymp(maxres) 

real*8 fc_poten(maxres) 

real*8 p_frag 

real*8 p^cap 

real*8 p_sul 

integer i 
integer ires 



get the coordinates of small molecule 

cal 1 readl i gand (atom_symbol , 1 i gand , 1 i gand_atomnum , 
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1 i gand_charge , maxsi ze) 



c* construct Z matrix of the interaction of full system *c 

c call exaGauss (coord, res_atomname,numres, res_atomnum, 

c $ maxres , mxratm , -1 , atom_symbol , 1 i gand , 

c $ 1 i gand-atomnum , maxsi ze , 500 , Theory) 



c* cut proteins into N pieces and N-1 caps *c 

cal 1 cut_protei n (coord , res^atomname , res_atomnum , res^name , 
$ numres , Frag , Frag^atomname , Frag_atomnum , 

$ cap , Cap_atomname , Cap_atomnum , maxres , 

$ mxratm , 1 evel , cut) 



c*--determine the #s of cys(CYX) residues & of di sulfur bond--*c 



cal 1 di sul f_bond (coord , res_atomname , res_atomnum , res_name , 
$ numres , Frag , Frag_atomname , Frag_atomnum , 

$ Sul cap , Sul cap_atomname , Sul Cap_atomnum , 

$ maxres , mx ratm , sn) 

c print*, '*********• ,sn, ' disulfur bond(s) found ********** 

c print* 



c* determine the charge for each piece and cap *c 

cal 1 get_pccharge (charge , endcharge , pcharge , ccharge , 
$ res„name , numres , maxres , 1 evel , cut) 

c* determine the polarity of fragment & cap *c 

cal 1 pol ari ty (res^name , Frag_i d , Cap_i d , numres , 
$ maxres , 1 evel , cut) 

c* calculate ab initio energies for each fragment *c 

cal 1 getenergy (Frag , Frag_atomname , Frag_atomnum , Fraa i d , 
$ pcharge, maxres, mxratm, atom^symbol ,ligand, 

$ 1 i gana.atomnum , 1 i gand_charqe , maxsi ze , 

$ numres , Theory , 1 , f_asymp , p_f rag , f gra_l i g) 

c* calculate ab initio energies for each caps *c 

cal 1 getenergy (cap , cap_atomname , Cap_atomnum , cap_i d , 
$ ccharge , maxres , mxratm , atom_symbol , 1 i gand , 

$ 1 i gana_atomnum , 1 i gand_charge , maxsi ze , 

$ numres , Theory , 2 , c_asymp , p_cap , eg ra_l i g) 

c* calculate ab initio energies for disulfur cap(s) *c 

if(sn.ge.l) then 
do 1=1, sn 

scharge(i)=0 
s_asymp(i)=0.dO 
end do 

cal 1 getenergy (Sul cap , Sul cap_atomname , sul Cap„atomnum , 
$ -l,scharge, maxres, mxratm, atom_symbol , 

$ 1 i gand , 1 1 gand_atomnum , 1 i gand_charge , 

$ maxsi ze,sn, Theory, l,s_asymp,p_sul ,sgra_lig) 

end if 



c* calculate total interaction energy *c 
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mam program.txt 

•^■^ 'M MFFFFF CCC 

■ ■^-'MMMMF C C 

" 'MMMMFFFF C C 

' 'M M M F C C 

■ 'M M F CCC 

*• "Normalized interaction energy == ' 

-p.cap-p_sul ) *27 . 2114d0*23 . 0605dO , 
print*, 'End of calculation' 



print*, 
print*, 
print*, 
print*, 
print*, 
print*, 



CCC 



CCC 
(p_fraq 



kcal/mol ' 



return 
end 



subrouti ne getenergy (tmp , tmp_atomname , tmp_atomnum , tmp_i d , 
$ charge , maxres , mxratm , atom_symbol , 1 i gand , 

$ 1 i gand_atomnum , 1 i gand_charge , maxsi ze , 

$ numres , Theory , start , asymp , abc , gradi ent) 

integer maxres 
integer mxratm 

real*8 tmp(maxres , mxratm, 3) 

character (len=4) tmp_atomname (maxres, mxratm) 

integer tmp_atomnumCmaxres) 

integer tmp_id (maxres) 

integer maxsi ze 

character (len=4) atom_symbol (maxsi ze) 
real*8 1 i gand (maxsi ze, 3) 
integer ligand_atomnum 
integer ligand^charge 

character (len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

real*8 gradi ent (maxsi ze, 3) 
integer charge(maxres) 
real*8 asympCmaxres) 
real*8 pes(maxres) 
real*8 abc, p_ligand 

integer id 
integer start 
real*8 min 
real*8 energy 
real*8 asyme 



calculate isolated ligand ab initio energy *c 

call get_unit(iout) 

cal 1 Gaussi an (tmp , tmp_atomname , tmp_atomnum , 1 , i out , 
$ charge (1) , max res, mxratm, atom^symbol , 

$ 1 i gand , 1 i gand_atomnum , 1 i gand_cnarge , 

$ maxsi ze , 2 , Theory) 

cal 1 Abi ni cal (i out , p_l i gand , gradi ent , 0 , maxsi ze , number*3 , 
$ (number+l i gand_atomnum) *3 , Theory) 
calculate the normalized interaction energy *c 

do ires=5tart, numres 

cal 1 mi ndi s (tmp , tmp_atomnum , tmp^atomname , i res , numres , 
$ 1 i gand , 1 i gand_atomnum , maxsi ze , maxres , 

$ mxratm, min) 
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cal 1 sel ectg roup (tmp_i d , i res , 4000 . 0 , 4200 . 0 , 4600 . 0 , 
$ min,maxres,id) 
if(id.ne.O) then 

call get_unit(iout) 

cal 1 Gaussi an (tmp , tmp_atomname , tmp_atomnum , i res , 
$ iout,cnarge(ires) ,fnaxres,mxratm,atonLsymbol , 

$ 1 i gand , 1 i gand_atomnum , 1 i gand_charge , maxsi ze , 

$ 1, Theory) 

number=tmp_atomnum(i res) 

cal 1 Abi ni cal (i out , energy , gradi ent , 0 , maxsi ze , number*3 , 
$ (number+1igand_atomnum)*3, Theory) 

pes(i res)=energy 

* calc. the corresponding isolated frag/cap energy once — *c 

if (asympCi res) .eg.O.dO) then 
call get_unit(iout) 

cal 1 Gaussi an (tmp , tmp„atomname , tmp„atomnum , i res , 
$ iout, charge(i res) ,maxres,mxratm, atom-Symbol , 

$ 1 i gand , 1 i gand_atomnum , 1 i gand_charge , maxsi ze , 

$ 0, Theory) 

cal 1 Abi ni cal (i out , energy , gradi ent , 0 , maxsi ze , number- 
$ (number+l i gand_atomnum) *3 , Theory) 

end if 

pes(i res)=pesCi res) -energy~p_li gand 
abc=abc+pes(i res) 
writeC*,*) ires, pes(ires) 
end if 
end do 

return 
end 

subrouti ne Abi ni cal (i out , energy , gradi ent , i d , maxsi ze , start , 
$ end, Theory) 

integer maxsi ze 
real*8 energy 

real*8 gradient(maxsize,3) 
integer id 
integer start 
integer end 
integer iout 
integer input^unit 
character(len=50) method 
character (len=100) grep 
character (len=100) cutl 
character (len=100) cut2 
character (len=3) tmp 
character (len=40) Theory 
integer i 
integer k 
integer kk 
integer kkk 

jout=iout 

do i=l. 3 

kkk=jout/10**2 

kk=(jout-kkk*100)/10 

k=(jout-kkk*100-kk*10) 
end do 

tmp=char(48+kkk)//char(48+kk)//char(48+k) 

Page 6 



main program.txt 
method='g98<fort. V/char(48+iout)//'> out V/tmp 

call System(method) 

if(Theory(l:3).ne. 'MPZ') then 

grep='grep "SCF Done" out V/tmp//'>100 V/tmp 
call System(grep) 

cutl='cut -d"=" -f2 100 V/tmp// '>200 V/tmp 
call System (cut 1) 

cut2='cut -d"A" -fl 200 V/tmp// VlOO V/tmp 
call System(cut2) 

call get_unit(input_unit) 

open (i nput_uni t , f i 1 e= * 100 ' //tmp , status= ' ol d ' ) 
read(input_unit,*) energy 
close(input_unit) 

else 

grep='grep "EUMP2" out V/tmp// •>100 V/tmp 
call System(grep) 

cutl='cut -d"-" -f3 100 V/tmp// •>200 V/tmp 
call System (cut 1) 

call get_unit(input_unit) 

open (i nput_uni t , f i 1 e=: ' 200 ' //tmp , status= ' ol d ' ) 
read(input_unit,*) energy 
cl ose (i nput_uni t) 

end if 

if(id.eq.l) then 

cal 1 f i ndgradi ent (gradi ent , start , end , maxsi ze) 
end if 

return 
end 



subroutine er rcorr (coord, res^atomname, res_atomnum, charge, 

endcharge , maxres , mxratm , atom_symbol , 1 i gand , 
1 i gand_atomnum , 1 i gand.charge , maxsi ze , 
Theory, level ,error,fc_poten) 

integer maxres ^ 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

character (1 en=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Cor r (maxres, mxratm, 3) 

cha ract e r (1 en=4) Co r r_at omname (max res , mx r atm) 

integer Corr_atomnum(maxres) 

integer Corr_charge(maxres) 

integer maxsi ze 



character (len=4) atom_symbol (maxsi ze) 
real*8 ligand(maxsize,3) 
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integer ligand_atomnum 
integer ligand_charge 

character(len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

integer charge(maxres) 
integer endcharge(2) 

integer ires 

integer iout 

real*8 fc_poten(maxres) 

real*8 energy 

do ires=l, numres-1 

call get_unit(iout) 

call combi neunits (Coord, res_atomname, res_atomnum, 
$ cor r , Cor r_atomname , Cor r_atomnum , i res , 

$ 1 evel , charge , endcharge , Cor recharge , 

$ maxres , mxratm , numres; 

cal 1 Gau s si an (Co r r , Co r r_a t omname , Co r r_atomnum , i r es , 
$ iout, Cor r_charge(i res) , maxres, mxratm, 

$ atom^symbol , 1 i gand , 1 i gand^atomnum , 

$ 1 i gand_charge , maxsi ze , 1 , Theory) 

call Abini cal (energy) 

f c_poten(i res)=energy 
end do 

return 
end 



subrouti ne f i ndgradi ent (gradi ent , start , end , maxsi ze) 

integer maxsi ze 
integer start 
integer end 
integer input_unit 
character ( I en=128) string 
character (len=20) var 
real*8 old 
real*8 der 

real*8 gradi ent (maxsi ze, 3) 
integer n 
logical s_eqi 

call get_unit(input_unit) 

open (uni t=i nput_uni t , f i 1 e= ' out ' , status= ' ol d ' ) 
do 

read(input_unit, ' (a) ') string 
if(s.eqi(string(2:9), 'Variable')) then 

read(input_unit, ' (a) ') string 

read(input_unit, ' (a) ') string 

do i=l, start 

read(input_unit, ' (a) ') string 

end do 

do i=start+l, end 
n=n+l 

read (St ri ng , ' (all , f 10 . 5 , f 10 . 5) ' ) var , ol d , der 
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gradient(l+(n-l)/3,n-3*C(n-l)/3))=der 
read(input_unit , ' (a) ■) string 
end do 
return 
end if 
end do 

cl ose (i nput^uni t) 

return 
end 

subroutine cut_protein(coord, res_atomname , res_atomnum, 
$ res_name , numres , Frag , Frag_atomname , 

$ Frag_atomnum , Cap , Cap_atomname , 

$ Cap_atomnum,maxres,mxratm, level ,cut) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

character (len=4) res^name (maxres) 

integer res.atomnum(maxres) 



real*8 Frag(maxres, mxratm, 3) 

character (len=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnumCmaxres) 

integer Frag_id(maxres) 

real*8 CapCmaxres, mxratm, 3) 
character(len=4) cap^atomname (maxres, mxratm) 
integer Cap_atomnum(maxres) 
integer Cap_id(maxres) 

character(len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

integer i res 

do ires=l, numres 

c* for the first residue in the protein chain *c 

if(ires.eq.l) then 

cal 1 resnol(coord , res_atomname , res.atomnum , 
$ Frag , Frag_atomname, Frag^atomnum, 

$ res.name , i res , max res , mxratm , 1 evel , cut) 



c* for the last residues in the protein chain *c 

else if (ires. eq. numres) then 

cal 1 resend(coord , res_atomname , res_atomnum , 
$ Frag , Frag^atomname , Frag_atomnum , 

$ res_name, ires, maxres, mxratm, level ,cut) 



c* for the middle residues in the protein chain- — 

else 

call resmi d (coord, res^atomname, res_atomnum, 
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$ Frag, Frag_atomname, Frag^atomnum, res_name, 

$ i res , numres , maxres , mxratm , 1 evel , cut) 

end if 



c* locate the caps in the middle of the protein chain *c 

if (1. It. ires) then 

call findcapCcoord, res_atomname, res_atomnum, 
$ Cap , Cap_atomname , Cap_atomnum , res_name , 

$ i res , numres , maxres , mxratm , 1 evel , cut) 

end if 
end do 

return 
end 

subrouti ne resnol(coord , res_atomname , res_atomnum , 
$ Frag , Frag_atomname , Frag_atomnum , 

$ res_name,i res, maxres, mxratm, level ,cut) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

character (len=4) res_name (maxres) 

integer res_atomnum (maxres) 



real*8 Frag (maxres, mxratm, 3) 

character ( I en=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnum(maxres) 

integer Frag_id (maxres) 

character (len=40) Theory 

integer level 

integer error 

integer cut 

integer numres 

integer ncpu 

integer i 

integer j 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

realms y(3) 

real*8 yl(3) 

character (len=4) Al 

character (len=4) A2 

character (len=4) AA 



copy the current residue: center of fragment 

call current_residue(coord, res_atomname, res_atomnum, 
$ Frag , Frag_atomname , Frag^atomnum , 

$ i res , maxres , mxratm) 
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c* construct the end cap(s) *c 

step=0 

cssssssssssssssssss cut=0 i.e. cut CA-N bond ssssssssssssssssssc 
if(cut.eq.O) then 

if (res_name(i res+1) .ne. 'PRO') then 
if (level .eq.l) then 

cal 1 CHSnext (coord , res^atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 

i f (1 evel . eq . 2 . or . 1 evel . eq . 3) then 

call CH2Rnext (coord, res_atomname, res_atomnum, ires, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 

if (level .eq. 4. or. level .eq. 5) then 

call CONHnext (coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
else 

cal 1 CH2Rnext (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

cssssssssssssssssss cut=l i.e. cut CA-C bond ssssssssssssssssssc 
if (cut. eq.l) then 

if (res_name(i res+1) .ne. 'PRO') then 

cal 1 NH2next (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

cal 1 CH2Rnext (coord , res_atomname , res_atomnum, i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

cssssssssssssssssss cut=2 i.e. cut C-N bond sssssssssssssssssssc 
if(cut.eq.2) then 

if (res_name(i res+l) .ne. 'PRO') then 
i f (1 evel . eq . 1 . or . 1 evel . eq . 3) then 

cal 1 CH3next (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

cal 1 CH2Rnext (coord , res^atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ max res , mxratm , step) 

end if 
else 

cal 1 CH2Rnext (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

csssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 

return 
end 
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subroutine resend (coord, res^atomname, res_atomnum, 
$ Frag , Frag_atomname , Frag_atomnum , 

$ res^name , i res , maxres , mxratm , 1 evel , cut) 

integer maxres 
integer mxratm 

real*8 coord (max res, mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

character (len=4) res^name (maxres) 

integer res_atomnum(maxres) 



real*8 Fraq(maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

i nt ege r F rag_atomn urn (max res) 

integer Frag_id(maxres) 

character (len=40) Theory 

integer level 

integer error 

integer cut 

integer numres 

integer ncpu 

integer i 

integer j 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

real*8 y(3) 

real*8 yl(3) 

character (len=4) Al 

character (len=4) a2 

character (len=4) AA 



c* copy the current residue: center of fragment *c 

call current_residue(coord, res_atomname, res_atomnum, 
$ Frag, Frag_atomname, Frag_atomn urn, 

$ i res , maxres , mxratm) 



c* construct the end cap(s) *c 

step=0 

cssssssssssssssssss cut=0 i.e. cut CA-n bond ssssssssssssssssssc 
if(cut.eq.O) then 

if (level .eq.l. or. level .eq. 2. or. level .eq.4) then 

call NH2prev(coord, res^atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum, i res , 

$ maxres , mxratm , step) 

end if 

i f (1 evel . eq . 3 . or . 1 evel . eq . 5) then 

call NHCOprev(coord, res^atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ max r es , mx ratm , step) 

end if 
end if 

cssssssssssssssssss cut=l i.e. cut CA-C bond ssssssssssssssssssc 
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if(cut.eq.l) then 

if (res_name(i res-l) .ne. 'PRO') then 
ifdevel .eq.l) then 

cal 1 CH3prev (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ max r es , mx ratm , s tep) 

end if 

if (level .eq.2) then 

call CH2Rprev (coord, res^atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 

if (level .eq.3) then 

cal 1 CHRNH2prev(coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres, mxratm, step) 

end if 

if (level .eq. 4. or .level .eq.5) then 

call CHRNHCOHprev(coord, res_atomname, res_atomnum,ires, 
$ Frag , Frag^atomname , Frag^atomnum, i res , 

$ maxres, mxratm, step) 

end if 
else 

cal 1 CHRNHCOHprev(coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum, i res , 

$ maxres , mxratm , step) 

end if 
end if 

cssssssssssssssssss cut=2 i.e. cut c-N bond sssssssssssssssssssc 
if (cut. eq.2) then 

if (res_name(i res-l) .ne. 'PRO') then 
i f (1 evel . eq . 1 . or . 1 evel . eq . 2) then 

cal 1 CH3prev(coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

call CH2Rprev(coord, res__atomname, res_atomnum,i res, 
$ Frag , Frag^atomname, Frag_atomnum, i res , 

$ maxres , mxratm , step) 

end if 
else 

call CHRNHCOHprev(coord,res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

csssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 

return 
end 



subroutine resmi d (coord, res_atomname, res^atomnum, 
$ Frag , Frag_atomname , Frag_atomnum , res_name , 

$ i res, numres, maxres, mxratm, level ,cut) 

integer maxres 
integer mxratm 



real*8 coord(maxres, mxratm, 3) 
character (1 en=4) res_atomname (maxres , mxratm) 
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character(len=4) res_name(maxres) 
integer res_atomnum(maxres) 



real*8 Fraq(maxres,mxratm,3) 

cha ract e r ( I en=4) Frag_atomname (max res , mx ratm) 

integer Frag_atomnumtmaxres) 

integer Frag_id(maxres) 

character (len=40) Theory 

integer level 

integer error 

integer cut 

integer numres 

integer ncpu 

integer i 

integer j 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

real*8 y(3) 

real*8 yl(3) 

character (len=4) Al 

character (len=4) a2 

character (len=4) AA 



copy the current residue: center of fragment 

cal 1 current^resi due (coord , res^atomname , res_atomnum, 
$ Frag , Frag^atomname , Frag_atomnum , 

$ i res , max res , mxratm) 



construct the end cap(s) *c 

step=0 

cssssssssssssssssss cut=0 i.e. cut CA-N bond ssssssssssssssssssc 
if(cut.eq.O) then 

if (level .eq.l) then 

call NH2prev(coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag^atomname , Frag_atomnum , i res , 

$ maxres , mx ratm , step) 

if (res_name(i res+1) .ne. ' PRO') then 

call CH3next (coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

call CH2Rnext (coord, res_atomname, res_atomnum,ires, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ max res , mx ratm , step) 

end if 
end if 

ifdevel .eq.2) then 

call NH2prev(coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

cal 1 CH2Rnext (coord , res_atomname , res^atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum, i res , 

$ max res , mxratm , step) 
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end if 

if (level .eq.3) then 

cal 1 NHCOp rev (coord , res^atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum , i res , 

$ maxres , mxratm , step) 

call CH2Rnext (coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 

if (level .eq.4) then 

cal 1 NH2prev(coord , res^atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

if (i res.eq. (numres-1)) then 

call CHZRnext (coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag^atomnum , i res , 

$ maxres , mxratm , step) 

else 

call CONHnext (coord, res__atomname, res_atomnum,i res, 
$ Frag , Frag^atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

i f (1 evel . eq . 5) then 

cal 1 NHCOp rev (coord , res_atomname , res^atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

if (i res.eq. (numres-1)) then 

call CH2Rnext (coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum, i res , 

$ maxres , mxratm , step) 

else 

cal 1 CONHnext (coord , res_atomname , res_atomnum , i res , 
$ Frag, Frag_atomname, Frag_atomnum,i res , 

$ maxres , mxratm , step) 

end if 
end if 
end if 

cssssssssssssssssss cut=l i.e. cut CA-C bond ssssssssssssssssssc 
if(cut.eq.l) then 

if (res_name(i res+l) .ne. ' PRO') then 

call NH2next(coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

cal 1 CH2Rnext (coord , res_atomname , res_atomnum , i res , 
$ Frag, Frag_atomname, Frag_atomnum,i res, 

$ maxres , mxratm, step) 

end if 

if (res_name(i res-l) .ne. ' PRO') then 
if (level .eq.l) then 

cal 1 CH3prev(coord , res^atomname , res.atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum, i res , 

$ maxres , mxratm , step) 

end if 

if (level .eq.2) then 

cal 1 CH2Rprev(coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum , i res , 
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$ max res , mx r atm , step; 

end if 

ifdevel .eq.3) then 
if (ires. eq. 2) then 
cal 1 CH2Rp rev (coord , res_atomname , res_atomnum, i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

cal 1 CHRNH2prev(coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ max res , mx ratm , step) 

end if 
end if 

i f (1 evel . eq . 4 . or . 1 evel . eq . 5) then 
if (i res.eq.2) then 

cal 1 CH2Rprev(coord , res_atomname , res.atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

cal 1 CHRNHCOHprev(coord , res_atomname , res_atomnum, i res , 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres, mxratm, step) 

end if 
end if 
else if (i res.eq.2) then 

call CH3prev(coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag^atomnum , i res , 

$ maxres , mxratm , step) 

else if (i res.gt. 2) then 

cal 1 CHRNHCOHprev(coord , res_atomname , res^atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

cssssssssssssssssss cut=2 i.e. cut C-N bond sssssssssssssssssssc 
if(cut.eq-2) then 

if (res_name(ires4-l) .ne. 'PRO') then 
if (level -eq.l. or. level .eq.3) then 

cal 1 CH3next (coord , res^atomname , res_atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

else 

cal 1 CH2Rnext (coord , res_atomname , res^atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum , i res , 

$ maxres , mxratm , step) 

end if 
else 

call CHZRnext (coord, res^atomname, res_atomnum,ires, 
$ Frag , Frag_atomname , Frag_atomnum, i res , 

$ maxres , mxratm , step) 

end if 

if (res_name(i res-1) ,ne. 'PRO') then 
if (level .eq.l. or. level .eq.2) then 

cal 1 CH3prev (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum, i res , 

$ maxres , mxratm , step) 

else 

call CH2Rprev(coord, res^atomname, res_atomnum,ires, 
$ Frag , Frag_atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 
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end if 
else if (ires. eq. 2) then 

cal 1 CHSprev (coord , res_atomname , res^atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum, i res , 

$ maxres , mx ratm , step) 

else if (i res.gt.2) then 

cal 1 CHRNHCOHprev(coord , res.atomname , res^atomnum , i res , 
$ Frag , Frag^atomname , Frag_atomnum , i res , 

$ maxres , mxratm , step) 

end if 
end if 

csssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 

return 
end 

c=============================================================== 

subroutine findcap(coord, res_atomname, res_atomnum, 
$ Cap , Cap.atomname , Cap_atomnum , res^name , 

$ i res , numres , maxres , mxratm , 1 evel , cut) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (1 en=4) res^atomname (maxres , mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

character (len=4) Cap_atomname (maxres, mxratm) 

integer Cap_atomnum(maxres) 

integer Cap_id(maxres) 

character (len=40) Theory 

integer level 

integer error 

integer cut 

integer numres 

integer ncpu 

integer i 

integer j 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

real*8 y(3) 

real*8 yl(3) 

character(len=4) Al 

characterOen=4) a2 

character(len=4) AA 



step=0 

cssssssssssssssssss cut=0 i.e. cut CA-N bond ssssssssssssssssssc 
if(cut.eq.O) then 

ifdevel .eq.l) then 

if (res_name(i res) .ne. 'PRO') then 

call Cap_CH3 (coord, res_atomname, res_atomnum, 
$ i res , Cap , cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 
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else 

call Cap_CH2R(coord, res^atomname, res^atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

end if 

cal 1 Cap_NH2 (coord , res_atomnanie , res^atomnum , 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm, step) 

end if 

if (level .eq.2) then 

call Cap_CH2R(coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap^atomnum , 

$ max res , mx r atm , step) 

call Cap_NH2(coord, res^atomname, res^atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ max res , mxratm , step) 

end if 

if (level .eq. 3) then 

cal 1 Cap_CH2R(coord , res.atomname , res_atomnum , 
$ ires, Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

cal 1 Cap_NHCO(coord , res_atomname, res_atomnum, 
$ i res , Cap , Cap^atomname , Cap^atomnum , 

$ maxres, mxratm, step) 

end if 

if (level .eq. 4) then 

if (i res.eq.numres) then 

call cap_CH2R (coord, res^atomname, res_atomnum, 
$ ires, Cap , Cap_atomname , Cap_atomnum , 

$ max res , mx ratm , step) 

else 

call Cap_CONH(coord, res^atomname, res_atomnum, 
$ i res , Cap , cap_atomname , Cap_atomnum , 

$ maxres, mxratm, step) 

end if 

call Cap_NH2 (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap^atomnum , 

$ max res , mxratm , step) 

end if 

if (level .eq.5) then 

if (i res.eq.numres) then 

call Cap_CH2R(coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

else 

call Cap_CONH(coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres, mxratm, step) 

end if 

cal 1 Cap_NHCO(coord , res^atomname , res^atomnum , 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

end if 
end if 

cssssssssssssssssss cut=l i.e. cut CA-C bond ssssssssssssssssssc 
if(cut.eq.l) then 

if(res_name(i res) .ne. 'PRO') then 

cal 1 ProCap_NH2 (coord , res_atomname , res^atomnum , 
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$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ max res, mx rat m, step) 

else 

call Cap_CH2R(coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

end if 

if (res_name(ires-l) .ne, 'PRO') then 
if (level .eq.l) then 

call ProCap_CH3 (coord, res^atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap^atomnum , 

$ maxres , mx ratm , step) 

end if 

if (level .eq.2) then 

cal 1 ProCap_CH2R(coord , res.atomname , res^atomnum , 
$ i res , Cap , Cap_atomname , Cap.atomnum , 

$ maxres, mxratm, step) 

end if 

if (level .eq.3) then 
if (ires. eq.2) then 

call ProCap_CH2R(coord, res^atomname, res_atomnum, 
$ i res , Cap , Cap^atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

else 

cal 1 ProCap_CHRNH2 (coord , res_atomname , res^atomnum , 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres, mxratm, step) 

end if 
end if 

i f (1 evel . eq . 4 . or . 1 evel . eq . 5) then 
if (ires. eq.2) then 

cal 1 ProCap_CH2R(coord , res_atomname , res^atomnum , 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

else 

call ProCap_CHRNHCOH (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap^atomnum , 

$ maxres, mxratm, step) 

end if 
end if 
else if (i res. eq.2) then 

cal 1 Procap_CH3 (coord , res^atomname , res_atomnum , 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ CapC , CapC_atomname , CapC^atomnum , 

$ max res , mx ratm , step) 

else if(ires.gt.2) then 

call ProCap_CHRNHCOH (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

end if 
end if 

cssssssssssssssssss cut=2 i.e. cut c-N bond sssssssssssssssssssc 
if(cut.eq-2) then 

i f ( res^name (i res) . ne . ' pro ' ) then 
ifdevel .eq.l. or .level .eq.3) then 

call Cap_CH3 (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap^atomname , Cap_atomnum , 

$ max res , mxratm , step) 

else 
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call Cap_CH2R(coord, res_atomname, res^atomnum, 
$ i res , cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

end if 
else 

call Cap_CH2R(coord, res_atomname, res.atomnum, 
$ i res , Cap , cap^atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

end if 

if (res_name(i res-1) .ne. ' PRO') then 
i f (1 evel . eq . 1 . or . 1 evel . eq . 2) then 

call ProCap_CH3 (coord, res_atomname, res_atomnum, 
$ i res , Cap , cap_atomname , Cap_atomnum , 

$ max r es , mx ratm , s t ep) 

else 

cal 1 ProCap_CH2R(coord , res_atomname , res_atomnum , 
$ i res , Cap , cap_atomname , Cap_atomnum , 

$ max res , mxratm , step) 

end if 
else if (ires. eq. 2) then 

call ProCap_CH3 (coord, res_atomname, res_atomnum, 
$ ires, Cap , cap_atomname , Cap^atomnum , 

$ max re s , mx ratm , step) 

else if(ires.gt.2) then 

cal 1 ProCap_CHRNHCOH(coord , res.atomname , res_atomnum , 
$ i res , Cap , Cap.atomname , Cap_atomnum , 

$ max res , mxratm , step) 

end if 
end if 

csssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 
Cap_atomnum(i res)=step 

return 
end 

0=============================================================== 

subrouti ne di sul f_bond (coord , res_atomname , res_atomnum, 
$ res_name,numres, Frag, Frag^atomname, Frag_atomnum, 

$ Sul Cap , sul Cap_atomname , Sul Cap^atomnum , 

$ maxres , mxratm , sn) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

character (len=4) res^name (maxres) 

integer res_atomnum(maxres) 



real*8 Fraq(maxres, mxratm, 3) 

character (Ten=4) Frag^atomname (max res, mxratm) 

integer Frag_atomnumCmaxres) 

integer Frag_id(maxres) 

real*8 SulCap(maxres, mxratm, 3) 

character (len=4) Sul cap_atomname(maxres, mxratm) 

integer SulCap_atomnum(maxres) 

integer sn 

integer numres 
integer num 
real*8 xl(3) 
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real*8 x2(3) 
real*8 dis 

real*8 sulfur(maxres,3) 
integer flag(maxres) 

sn=0 
num=0 

do ires=l, numres 

if (res_name(i res) .eq. 'CYX' .or. res_name(i res) .eq. 'CYS') then 
do iatom=l, res_atomnum(i res) 

if (res_atomname(ires,iatom).eq. 'SG') then 
num=num+l 
flaq(num)=ires 
do k=l, 3 

sul f ur (num, k)=coord(i res , i atom, k) 
end do 
end if 
end do 
end if 
end do 

if (num. ge. 2) then 
do i=l, num 
do j=l, num 

if (i .ne. j .and.i .It . j) then 
do k=l, 3 

xl(k)=sulfur(i ,k) 
x2(k)=sulfur(j,k) 
end do 
dis=0.dO 

cal 1 di stanceCxl , x2 , di s) 

if (dis.lt.2.5dO) then 

call coordupdate (coord, res_atomname, res_atomnum, 
$ Frag , Frag_atomname , Frag_atomnum , 

$ Sul Cap , Sul Cap^atomname , Sul Cap_atomnum , 

$ flag(i) ,flag(j) ,maxres,mxratm,sn) 

end if 
end if 
end do 
end do 
end if 

return 
end 

subroutine coordupdate (coord, res_atomname, res_atomnum, 
$ Frag , Frag_atomname , Frag_atomnum , 

$ Sul Cap , Sul Cap_atomname , Sul Cap^atomnum , 

$ m , n , maxres , mxratm , sn) 

integer maxres 
integer mxratm 

real*8 coord(maxres , mxratm, 3) 

character (len=4) res^atomname (maxres, mxratm) 

character (len=4) res_name (maxres) 

integer res_atomnum(maxres) 



real*8 Frag (maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnumCmaxres) 
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integer Frag_id(maxres) 

real*8 SulCap(maxres,mxratm,3) 

character (1 en=4) Sul Cap_atomname(maxres , mxratm) 

integer SulCap_atomnum(maxres) 

integer sn 

integer i 

integer q 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

real*8 y(3) 

real*8 yl(3) 

character(len=4) Al 

character (len=4) a2 

character (len=4) AA 



real*8 tmpl(maxres, mxratm, 3) 
character (len=4) tmpl_atomname(maxres, mxratm) 
rea1*8 tmp2(maxres, mxratm, 3) 
character(len=4) tmp2_atomname(maxres, mxratm) 

integer m 
integer n 

integer flag(maxres) 

do i=l, 2 

if(i.eq.l) then 

i res=m 

jres=n 
else 

i res=n 

j res=m 
end if 

Al='CB' 
A2='CA' 
AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
$ j res , j res , Al , a2 , maxres , mxratm , x , y) 

call bondlength(x,y,AA) 
step=0 

do iatom=l, res_atomnum(jres) 

if (res_atomname( jres, 1 atom) .eq. 'SG' .or . 
$ res_atomname ( j res , i atom) . eq . ' CB ' . or . 

$ res_atomname ( j res , i atom) . eq . * HB " . or . 

$ res_atomname (] res , i atom) . eq . * CA ' ) then 

step=step+l 
do I<=1, 3 

i f (i . eq . 1) tmpl(i res , step , k)=coord( j res , i atom , k) 
i f (i . eq . 2) tmp2 (i res , step , k) =coord ( j res , i atom , k) 
end do 

if(res_atomname(j res, i atom) -eq. 'CA') then 
i f (i . eq . 1) tmpLatomname(i res , step)= " H ' 
i f (i . eq . 2) tmp2_atomname(i res , step)= " H ' 
do k=l, 3 

if(i .eq.l) tmpl(i res, step, k)=y(k) 
if(i .eq.2) tmp2(i res, step, k)=y(k) 
end do 
else 
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i f (i . eq . 1) tmpl^atomname (i res , step)= 
$ res^atomname ( j res , i atom) 

i f (i . eq . 2) tmp2_atomnameO res , step)= 
$ res_atomname ( j res , i atom) 

end if 
end if 
end do 
end do 

do i=l, step 
do k=l, 3 

Frag(m,Frag_atomnum(m)+i ,k)=tmpl(m,i ,k) 
Frag(n,Frag_atomnum(n)+i ,k)=tmp2(n,i ,k) 
end do 

Frag_atomname(m, Frag_atomnum(m)+i )=tmpLatomname(m, i ) 
Frag_atomname(n , Frag_atomnum(n)+i )=tmp2_atomname(n , i ) 
end do 

Frag_atomnum(m)=Frag_atomnum(m)+step 
Frag_atomnum(n)=Frag_atomnum(n)+step 

sn=sn+l 

do i=l, step 

do k=l, 3 

SulCap(sn,i ,k)=tmpl(m,i ,k) 

end do 

sul Cap^atomnameCsn , i )=tmpl_atomname(m , i ) 
end do 

do i=l, step 
do k=l, 3 

Su1cap(sn,step+i ,k)=tmp2(n,i ,k) 
end do 

Sul Cap_atomname (sn , step+i )=tmp2_atomname (n , i ) 
end do 

SulCap_atomnum(sn)=step+step 

return 
end 



subroutine readligand(atom_symbol ,ligand,ligand«atomnum, 
$ In gand_charge , maxsi ze) 

integer maxsi ze 

character (1 en=4) atom_symbol (maxsi ze) 
real*8 ligand(maxsize,3) 
integer ligand_atomnum 
integer ligand_charge 



open (111 , f i 1 e= ' 1 i gand . dat ' , status= ' ol d ' ) 
read (111,*) ligand_atomnum, ligand_charge 
do i=l, 1igand_atomnum 

read (111 , *) atom_symbol (i ) . (1 i gand (i , j ) , j =1 . 3) 
end do 
close(lll) 

return 
end 



subrouti ne exaGauss (coord , res_atomname , numres , res_atomnum, 
$ maxres , mxratm , charge , atom_symbol , 1 i gand , 

$ 1 i gand_atomnum , maxsi ze , i out , Theory) 
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integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (len=4) res.atomname (maxres, mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



integer maxsize 

character (len=4) atom_symbol (maxsize) 
real*8 ligand (maxsize, 3) 
integer 1 i gand_atomnum 
integer ligand^charge 

character (len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

integer charge 
integer ires 
integer iatom 
integer iout 

wri te (i out , " ( ' #T ' , 2x , a20 , 2x , ' SCF(MaxCycl e=500 , 
$ converts) test')") Theory 

write(iout , *) 

write(iout, "Cab initio calculation for full system')") 
write(iout,*) 

write(iout, "(12, Ix, •!•)") charge 

do ires=l, numres 

do iatom=l, res_atomnum(i res) 

wri te (i out , 111) res^atomname (i res , i atom) (1:1), 
$ coord (i res , i atom ,1:3) 

111 format (a , 2x , 3gl4 . 6) 

end do 
end do 

do i=l, ligand_atomnum 

write(iout,lll) atom^symbol (i ) , ligand(i ,1:3) 
end do 

write(iout,*) 
close(iout) 

return 
end 



subrouti ne Gaussi an (tmp , tmp^atomname , tmp_atomnum , i res , i out , 
$ charge , maxres , mxratm , atom_symbol , 1 i gand , 

$ 1 i gand_atomnum , 1 i gand^charge , maxsi ze , i d , 

$ Theory) 

integer maxres 
integer mxratm 

real*8 tmp(maxres, mxratm, 3) 
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character (len=4) tmp_atomname(maxres,mxratm) 
integer tmp_atomnum(maxres) 
integer tmp_id(maxres) 
integer maxsize 

character(len=4) atom_symbol (maxsize) 
real*8 1 i gand (maxsize, 3 j 
integer ligand_atomnum 
integer ligand_charge 

character (len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

integer tmp_charge 
integer id 
integer charge 
integer iout 
integer j 

tmp_charge=l i gand_charge 
numatm=tmp_atomnum(i res) 

write(iout,"('#T' ,2x,a21,2x, 
$ 'NoSymm SCF(MaxCycle=500,Conver=5) test')") Theory 

write(iout,*) 

wri te (i out , " ( ' 1 i gand+pepti de potenti al ' ) ") 
write(iout,*) 

if(id.ne.2) then 

if(id.eq.O) tmp_charge=0 

wri te (i out , " (I2 , Ix , ' 1 ' ) ") charge+tmp_charge 

do j=l, numatm 

wr i t e (i ou t , 111) tmp_atomname (i res , j ) , tmp (i res , j , 1 : 3) 
111 format (a4 , 2x , 3gl4 . 6) 

end do 

if(id.eq.l) then 

do 3=1, ligand^atomnum 

wri te (iout, 111) atom_symbol ( j ) , 1 i gand ( j , 1 : 3) 
end do 
end if 
else 

write(iout,"(l2,lx, '1')") tmp^charge 
do j=l, ligand^atomnum 

wri te (i out , 111) atom_symbol ( j ) , 1 i gand ( j , 1 : 3) 
end do 
end if 

wri te(iout , *) 
c write(iout,"('— linkl— ')") 

close(iout) 

return 
end 

c 

C==============================================:================= 

subrouti ne NH2next (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag^atomname , Frag^atomnum , j res , 
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max res , mxratm , step) 



integer maxres 
integer mxratm 

real*8 coord(maxres , mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

character Clen=4) res^name (maxres) 

integer res_atomnum(maxres) 



real*8 Frag (maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnum(maxres) 

integer Frag_id (maxres) 

integer i 
integer q 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character(len=4) Al 
character(len=4) A2 
character (len=4) AA 



A1='N' 
A2='CA' 
AA= • CN • 

cal 1 twopoi nts (coord , res_atomnum , res^atomname , 
$ i res+1 , i res+1 , Al , a2 , max res , mxratm , x , y) 

call bondlength(x,y,AA) 
step=0 

do iatom=l, res_atomnum(i res+1) 

if (res_atomname(i res+1, i atom) .eq. 'N' .or. 
$ res_atomname(i res+1, i atom) .eq. 'H' .or. 

$ res_atomname (i res+1, i atom) .eq. 'CA') then 

step=step+l 
do k=l, 3 

Frag(j res , Frag_atomnum(j res)+step , k)= 
$ coord (i res+1, i atom, k) 

end do 

i f(res_atomname(i res+1, i atom) .eq. 'ca') then 

Fraq_atomname( j res , Frag_atomnum( j res)+step)= ' H ' 
do k=l, 3 

Frag ( j res , Frag_atomnum( j res)+step , k)=y (k) 
end do 
else 

Frag_atomname(j res , Frag_atomnum(j res)+step)= 
$ res_atomname (i res+1 , i atom) (1 : 1) 

end if 
end if 
end do 

Frag_atomnum( j res)=Frag_atomnum( j res)+step 

return 
end 
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c============================================================== 

subroutine CH3prev(coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag.atomnum , j res , 

$ maxres,mxratm,step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

cha racte r (1 en=4) r es_atomname (max res , mx ratm) 

characterO en=4) res_hame(maxres) 

integer res_atomnum(maxres) 



rea1*8 Frag (maxres , mxratm, 3) 

character ( I en=4) Frag_atomname (maxres, mxratm) 

integer Frag_atomnum(maxres) 

integer Frag_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
characterO en=4) Al 
characterO en=4) a2 
characterO en=4) aa 



$ 
$ 
$ 
$ 
$ 
$ 
$ 
$ 



Al='CA' 

A2='N' 

AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
i res- 1 , i res-1 , Al , A2 , max res , mx ratm , x , y) 
call bondlength(x,y,AA) 

Al='CA' 
A2='CB' 

cal 1 twopoi nts (coord , res_atomnum , res^atomname , 

i res-1 , i res-1 , Al , a2 , maxres , mxratm , x , yl) 
cal 1 bondl ength (x , yl , AA) 
step=0 

do iatom=l, res_atomnum(i res-1) 

if (res_atomname(i res-1, i atom) .eq. 'c 



res_atomname(i i 
res_atomname(i i 
res_atomname(i 
res_atomname(i 
res_atomname(i 
res_atomname(i 
res_atomname(i 



'-or. 

i res-1, iatbm) .eq. '0' .or. 
i res-1, i atom) .eq. 'CA' .or. 
i res-1, iatom) .eq. 'HA' .or. 
i res-1, i atom). eq. 'haI' .or. 
i res-1, iatom) .eq. 'CB' .or. 
1 res-1, iatom) .eq. •ha2' .or. 
i res-1, iatom) .eq. 'HA3' .or. 
res_atomname(i res-1, iatom). eq. 'N') then 
step=step+l 
do k=l, 3 

Frag( j res , Frag_atomnum( j res)+step , k)= 
coord(i res-1, i atom, k) 

end do 

if(res_atomname(i res-1, iatom) .eq. 'N') then 

Frag_atomname( j res , Frag_atomnum( j res)+step)= ' H ' 
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do k=l, 3 

Frag C j res , Frag_atomnum( j res)+step , k)=y (k) 
end do 
end if 

if (res_atomname(i res-l,iatom) .eq. 'CB') then 
FraaatomnameCj res , Frag_atomnum( j res)+step)= ' H " 
do k=l, 3 

Frag( j res , Frag_atomnum(j res)+step , k)=yl(k) 
end do 
end if 

if (res_atomname(i res-l,iatom) .ne. 'n' .and. 
$ res_atomname(ires-l,iatoni) .ne. 'CB') then 

Frag_atomname(j res , Frag_atomnum( j res)+step)= 
$ res_atomname(i res-1 , i atom) (1 : 1) 

end if 
end if 
end do 

Frag_atomnum( j res)=Frag_atomnum( j res)+step 

return 
end 

subrouti ne CH2Rprev (coord , res^atomname , res_atomnum, i res , 
$ Frag , Frag^atomname , Frag_atomnum , j res , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (1 en=4) res_atomname (maxres , mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



rea1*8 Fraq (maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnum(maxres) 

integer Frag^id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) AA 



A1='CA' 

A2='N' 

AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
$ i res-1 , i res-1 , Al , a2 , max res , mxratm , x , y) 

call bondlength(x,y,AA) 
step=0 
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do iatom=l, res_atomnum(ires-l) 

if (res_atoinname(i res-l,iatom) .ne. 'H' .and. 
$ res_atomname(ires-l,iatom) .ne. 'Hi' .and. 

$ res_atomnameCi res-l,iatom) .ne. 'H2' .and. 

$ res_atomname(i res-l,iatom) .ne. 'h3') then 

step=step+l 
do k=l, 3 

FragCj res , Frag_atomnum( j res)+step , k)= 
$ coord(i res-l,iatom,k) 

end do 

if (res_atomname(i res-l,iatom) .eq. 'N') then 

FraQ_atomname( j res , Frag_atomnum(j res)+step)= ' H ' 
do k=l, 3 

FragCj res , Frag_atomnum( j res)+step , k)=y (k) 
end do 
else 

Frag_atomname(j res , Frag_atomnum(j res)+step)= 
$ res_atomnanie(i res-1, i atom) (1:1) 

end if 
end if 
end do 

Frag_atomnum( j res)=Frag_atomnum( j res)+step 

return 
end 

subroutine CHRNH2prev(coord, res_atomname, res_atomnum, i res, 
$ Frag , Frag^atomname , Frag_atomnum , j res , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

rea1*8 coord(maxres, mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Frag (maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnum(maxres) 

integer Frag_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) a2 
character (len=4) aa 



Al='N' 
A2='C' 
AA= • CN • 
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cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
$ i res - 1 , i re s - 2 , Al , a2 , max r e s , mx r atm , x , y) 

call bondlength(x,y,AA) 
step=0 

do iatom=l, res_atomnum(i res-1) 
step=step+l 
do k=l, 3 

FragC j res , Frag_atomnum(g res)+step , k)= 
$ coord (i res-l , i atom , k) 

end do 

Frag_atomname(jres,Frag_atomnum(j res)+step)= 
$ res_atomname(i res-l, i atom) (1:1) 

end do 

do iatom=l, res_atomnum(i res-2) 

i f (res_atomname(i res-2 , i atom) . eq . ' C ) then 
step=step+l 
do k=l, 3 

Frag ( j res , Frag_atomnum( j res)+step , k)=y (k) 
end do 

Frag_atomname(j res , Frag_atomnum( j res)+step)= ' H ' 
end if 
end do 

Frag_atomnum(j res)=Frag_atomnum(j res)+step 

return 
end 



subroutine CHRNHCOHprev(coord, res_atomname, res_atomnum,i res , 
$ Frag , Frag_atomname , Frag^atomnum , j res , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

cha ract e r (1 en=4) r es_atomname (max r es , mx ratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Fraq(maxres, mxratm, 3) 

character (len=4) Frag__atomname (maxres, mxratm) 

integer Frag_atomnumCmaxres) 

integer Frag_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
realms y(3) 
realms yl(3) 
character(len=4) Al 
character (len=4) a2 
character (len=4) AA 



Al='c' 
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a2='<:a' 

AA='CC* 

cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
$ ire5-2,ires-2 , Al, A2 ,maxres ,mxratm, x,y) 

cal 1 bondl ength (x , y , AA) 
step=0 

do iatom=l, res_atomnum(i res-l) 
step=step+l 
do k=l, 3 

FragC j res , Frag_atomnum( j res)+step , k)= 
$ coord (i res- 1 , 1 atom , k) 

end do 

Frag_atomname(j res , Frag_atomnum(j res)+step)= 
$ res_atomname(i res-l, i atom) (1:1) 

end do 

do iatom=l, res„atomnum(i res-2) 

if (res_atomname(i res-2,iatom) .eq. 'C' .or . 
$ res_atomname(ires-2,iatom).eq. .or. 

$ res_atomname(ires-2,iatom) -eq. "CA') then 

step=step+l 
do k=l, 3 

Frag( j res , Frag_atomnum( j res)+step , k)= 
$ coo r d (i r es - 2 , i atom , k) 

end do 

if (res_atomname(i res-2,iatom) .eq. 'CA') then 

Frag_atomname(j res , Frag_atomnum( j res)+step) = ' H ' 
do k=l, 3 

Frag (j res, Frag_atomnum(jres)+step,k)=y(k) 
end do 
else 

Frag_atomname(j res , Frag_atomnum(i res)+step)= 
$ res_atomname(i res-2 , i atom) (1:1) 

end if 
end if 
end do 

Frag_atomnum(j res)=Frag_atomnum(j res)+step 

return 
end 

subroutine PrpCap_NH2 (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

cha racte r (1 en=4) r es^name (max res) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

char act e r (1 en=4) Cap_atomname (max res , mx ratm) 

integer Cap_atomnum(maxres) 

integer Cap_id (maxres) 

integer i 
integer j 
integer latom 
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integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) AA 



Al='N' 
A2='CA* 
AA= ' CN ' 

call twopoi nts (coord, res_atomnum, res_atomname , 

i res.i res,Al,A2,maxres,mxratm,x,y) 
call bondlength(x,y,AA) 

do iatom=l, res_atomnum(i res) 

if (res_atomname(i res,iatom) .eq. 'N' .or. 
res_atomname(i res,iatom) .eq. 'H* .or . 
res_atomname(i res, i atom) .eq. 'CA') then 
res_atomname(i res, i atom) .eq. 'CB') then 
step=step+l 
do k=l, 3 

Cap(i res , step , k)=coord (i res , i atom , k) 
end do 

if (res_atomname(ires,iatom) .eq. 'CA') then 

res_atomname (i res , i atom) . eq . ' CB ' ) then 

Cap_atomname(i res , step)= ' H ' 

do k=l, 3 

CapCi res , step , k)=y (k) 

end do 
else 

Cap_atomname(i res , step)= 

res_atomname (i res , i atom) (1 : 1) 

end if 
end if 
end do 

return 
end 



subroutine ProCap_CH3 (coord, res_atomname, res_atomnum, 
i res , Cap , cap_atomname , Cap^atomnum , 
maxres ,mxratm, step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

cha racte r (1 en=4) res_atomname (max res , mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Cap (maxres, mxratm, 3) 

characte r (1 en=4) cap_atomname (max res , mx ratm) 

integer Cap_atomnum(maxres) 

integer Cap_id (maxres) 

integer i 
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integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
rea1*8 x(3) 
rea1*8 y(3) 
real*8 yl(3) 
character(len=4) Al 
character (len=4) A2 
character (len=4) AA 



mam program.txt 



Al='CA' 

A2='N' 

AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res.atomname , 
$ i re s - 1 , i res - 1 , Al , a2 , max re s , mx ratm , x , y ) 

call bondlength(x,y,AA) 

Al='CA' 
A2='CB' 

call twopoi nts (coo rd, res_atomnum,res_atomname, 
$ i res-l,i res-l,Al,A2,maxres,mxratm,x,yl) 

cal 1 bondl ength (x , yl , AA) 

do iatom=l, res_atomnum(i res-1) 

if (res_atomname(i res-l,iatom) .eq. 'C' .or. 
$ res_atomname(i res-l,iatom) .eq, 'O' .or. 

$ res_atomname(ires-l,iatom) .eq. 'CA' .or. 

$ res_atomname(i res-l,iatom) .eq. 'HA' .or . 

$ res_atomname(i res-l,iatom) .eq. 'HAl' .or. 

$ res_atomname(i res-l,iatom) .eq. 'CB' .or . 

$ res_atomname(i res-l,iatom) .eq. 'HA2' .or. 

$ res_atomname(i res-l,iatom) .eq. 'HA3' .or. 

$ res_atomname(i res- l,i atom) .eq. 'N') then 

step=step+l 

do k=l, 3 

Cap(i res , step, k)=coord(i res-1, iatom, k) 
end do 

i f(res_atomname(i res-1, i atom) .eq. 'N') then 

Cap_atomname (i res , step)= ' H ' 

do k=l, 3 

cap(i res, step, k)=y(k) 

end do 
end if 

if (res_atomname(i res-1, iatom) .eq. 'CB') then 

Cap^atomname (i res , step)= ' H ' 

do k=l, 3 

Cap(i res , step , k)=yl(k) 

end do 
end if 

if (res_atomname(i res-1, iatom) .ne. 'N' .and. 
$ res_atomname(i res-1, iatom) .ne. 'CB') then 

Cap_atomname(i res , step)= 
$ res_atomname (i res-1 , i atom) (1:1) 

end if 
end if 
end do 



return 
end 
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subroutine ProCap_CH2R(coord, res_atomname, res^atomnum, 
ires, Cap , Cap_atomname , Cap_atomnum , 
maxres , mxratm, step) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (1 en=4) res_atomname(maxres , mxratm) 

characterOen=4) res_name(maxres) 

integer res_atomnum(maxres) 



rea1*8 Cap(maxres, mxratm, 3) 

character (len=4) Cap_atomname(maxres, mxratm) 

integer Cap_atomnum(maxres) 

integer Cap_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer i res 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character Oen=4) Al 
character Oen=4) a2 
character Oen=4) AA 



Al='CA' 

A2='N' 

AA='CC' 

cal 1 twopoi nts (coord , res^atomnum , res_atomname , 
i res- 1 , i res-1 , Al , A2 , maxres , mxratm , x , y) 
call bondlength(x,y,AA) 

do iatom=l, res_atomnum(i res-1) 

if (res_atomname(i res-1, i atom) .ne. 'H' .and. 

res_atomname(i res-1, i atom) .ne. 'Hi' .and. 
res_atomname(i res-1, i atom) .ne. 'H2' .and. 
res_atomname(i res-1, i atom) .ne. •h3') then 
step=step+l 
do k=l, 3 

Cap(i res , step, k)=coord(i res-1, i atom, k) 
end do 

if(res_atomname(i res-1, i atom). eq. 'N') then 

Cap_atomname(i res , step)= ' H ' 

do k=l, 3 

cap(i res , step , k)=y(k) 

end do 
else 

Cap_atomname(i res , step)= 

res_atomname(i res-1, i atom) (1:1) 

end if 
end if 
end do 

return 
end 
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c=============================================================== 

subroutine ProCap_CHRNH2 (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

character (len=4) res^name (maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

characte r (1 en=4) Cap_atomname (maxres , mx ratm) 

integer cap_atomnum (maxres) 

integer Cap_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) AA 



A1='N' 
A2='C' 
AA='CN' 

call twopoi nts (coord, res_atomnum, res_atomname, 
$ i r e s - 1 , i r es - 2 , Al , a2 , max r es , mx ratm , x , y ) 

call bondlength(x,y,AA) 

do iatom=l, res_atomnum(i res-1) 
step=step+l 
do k=l, 3 

Cap(i res , step , k)=coord(i res-1, i atom , k) 
end do 

Cap_atomname(i res , step)= 
$ res_atomname(i res-1, i atom) (1:1) 

end do 

do iatom=l, res_atomnum(i res-2) 

if (res_atomname(i res-2, i atom) .eq. 'C) then 
step=step+l 
do k=l, 3 

cap(i res , step , k)=y (k) 
end do 

cap^atomname (i res , step)= ' H ' 
end if 
end do 

return 
end 

c====================================================== 
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subrouti ne ProCap_CHRNHCOH (coord , res.atomname , res_atomnum, 
i res , cap , Cap_atomname , cap_atomnum , 
maxres , mxratm , step) 

integer maxres 
integer mxratm 

rea1*8 coord (maxres, mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

character Oen=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres , mxratm, 3) 

character (1 en=4) Cap_atomname (maxres , mxratm) 

integer Cap_atomnum(maxres) 

integer Cap_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) AA 



Al='C' 

A2='CA' 

AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res^atomname , 
i res-2 , i res-2 , Al , a2 , max res , mxratm , x , y) 
call bondlength(x,y,AA) 

do iatom=l, res_atomnum(i res-1) 
step=step+l 
do k=l, 3 

Cap(i res , step, k)=coord(i res-1, i atom, k) 
end do 

Cap_atomname(i res , step)= 

res_atomname(i res-1, i atom) (1:1) 

end do 

do iatom=l, res_atomnum(i res-2) 

if (res_atomname(i res - 2, i atom) .eq. 'C .or. 
res_atomname(i res-2 ,i atom) .eq. 'o* .or . 
res_atomname(i res-2, i atom) .eq. 'CA') then 
step=step+l 
do k=l, 3 

CapO res , step, k)=coord(i res-2 ,i atom, k) 
end do 

if (res_atomname(i res-2, i atom) .eq. 'CA') then 

Cap_atomname(i res , step)= ' H ' 

do k=l, 3 

Cap(i res , step, k)=y (k) 

end do 
else 
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Cap_atomname(i res , step)= 
$ res_atoinname(i res-2 , i atom) (1 : 1) 

end if 
end if 
end do 

return 
end 



subroutine readpdb(coord, res_atomname, res_atomnum, res_name, 
$ charge , endcharge , numres , maxres , mxratm) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

character (len=4) res_name (maxres) 

integer res_atomnum(maxres) 



integer charge(maxres) 

integer endcFiarge(2) 

integer atomnumCmaxres, mxratm) 

character (len=80) filepdb 



integer 
integer 
integer 
integer 
integer 
integer 



arg 
argc 
error 
len 

nput_uni t 

OS 



integer ipxfargc 

integer lenc 

integer numarg 

integer numres 



Open the pdb file *c 

call get_unit(input_unit) 

open (uni t=i nput_uni t , f i 1 e= ' protei n . pdb ' , status= ' ol d " , i ostat=i os) 
if (ios.ne.O) then 
writeC*,*) ' * 

writeC*,*) 'PDB_PRB - Fatal error!' 
writeC*,*) 'Could not open the PDB file.' 
stop 
end if 



* Read the information from the file *c 

call pdb_read (coord, res^atomname, res_atomnum, res_name, 
$ atomnum , i nput_uni t , numres , maxres , mxratm) 



close the pdb file *c 

cl ose (uni t=i nput_uni t) 



f-jncl charge for each residue *c 

call find_charge(res_atomname, res^atomnum, res_name, 
$ charge , endcharge , numres , maxres , mxratm) 
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c*--change coordinates according to the pdb file from chem3D--*c 
c cal 1 change_coord (coo rd , res_atomnum , atomnum , numres , 

c $ maxres.mxratm) 



c* find CMS *c 

c cal 1 f i nd^cms (coord , res_atomname , res^atomnum , numres , 

c $ maxres,mxratm) 



c* after CMS transformation, print out the pdb file *c 

c call print_pdb(coord, res_atomname, res_atomnum, res_name, 

c $ numres, maxres.mxratm) 

return 
end 

C==========================================================:===== 

subroutine get_unit(iunit) 

integer i 
integer ios 
integer iunit 
logical lopen 

iunit=0 

do i=l, 500 

if (i .ne. S.and.i .ne.6) then 

i nqui re (uni t=i , opened=l open , i ostat=i os) 
if(ios.eq.O) then 

if (.not .lopen) then 
iuni t=i 
return 
end if 
end if 
end if 
end do 

return 
end 

c=============================================================== 

subroutine pdb_init(coord, res_atomname, res_atomnum, 
$ numres , maxres , mxratm) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (len=4) res^atomname (max res, mxratm) 

character (len=4) res_name (maxres) 

integer res_atomnum(maxres) 



integer charge(maxres) 
integer numres 

coo rd (1 : max res , 1 : mx ratm , 1 : 3) =0 . 0 
res_atomnum(l:maxres)=0 
res_atomname (1 : maxres , 1 : mxratm)= ' ' 
charged: maxres)=0 
numres=0 

Page 38 



mam program.txt 

return 
end 

subroutine pdb_read (coord, res_atomname, res^atomnum, 

res_name , atomnum , i nput^uni t , numres , maxres , mxratm) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 
character(len=4) res_atomname(maxres, mxratm) 
character (len=4) res_name (maxres) 
integer res_atomnum(maxres) 



integer atomnum(maxres, mxratm) 

integer ratm 

integer ibase 

integer input_unit 

integer ios 

integer numres 

integer prvnumres 

character prvchn 

integer prvresno 

character (len=3) prvresname 

character (len=128) string 

integer iatom 
integer resno 
real*8 x 
real*8 y 
real*8 z 
real*8 occ 
real*8 temp 

character (len=4) wl 
character (len=4) atomname 
character altloc 
character (len=3) resname 
character chains 
character icode 
character (len=4) segid 
character (len=2) element 
character (len=2) charge 

logical s^eqi 



c* initialize PDB information 

call pdb_init(coord, res^atomname, res^atomnum, 
$ numres, maxres, mxratm) 

ibase=0 
prvchn=* ' 
prvresno=0 
prvresname=* ' 

do 

read(input_unit, ' (a) ' ,iostat=ios) string 
if(ios.ne-O) then 
exit 
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end if 

i f (s_eqi (st ri ng (1 : 6) , ' ENDMDL ' ) ) then 

exit 

else if(s_eqi(string(l:4), 'ATOM')) then 
readCstring, 



$ ' (a6 , i 5 , Ix , a4 , al , a3 , Ix , al , i 4 , al , 3x , 
$ 3f8.3,2f6.2,6x,a4,a2,a2)' , 

$ iostat=ios) 

$ wl , i atom , atomname , al 1 1 oc , resname , chai ns , resno , 

$ i code , X , y , z , occ , temp , segi d , el ement , charge 



if(ios.ne.O) then 

exit 
end if 



c* — Remove a possible initial blank in atomname or resname — *c 
if (atomname(l:l) .eq. ' ') then 

atomname=atomname(2 : ) 
end if 

if (resname(l:l) .eq. * ') then 

resname = resname(2:) 
end if 

if (atomname (1: 1) .eq. '1' .or. atomname (1:1) .eq. '2' .or. 
$ atomname(l:l) .eq. '3') then 

atomname=atomname(2 : ) 
end if 



c* If necessary, increment the number of residues read *c 

i f (resno . ne . prvresno . or . resname . ne . prvresname 
$ .or.chains.ne.prvchn) then 

prvresno=resno 
prvresname=resname 
numres=numres+l 
prvchn=chains 
end if 



c*--For each atom, store the atomic coordinate, and the name--*c 

c* and number of the residue to which the atom belongs *c 

i f (1 . 1 e . numres . and . numres . 1 e . maxres) then 
i f (numres . ne . prvnumres) ratm=0 

ratm=ratm-i-l 

res_name(numres)=resname 
res_atomnum(numres)=ratm 
coord (numres, ratm,l)=x 
coord (numres, ratm,2)=y 
coord (numres, ratm,3)=z 
res_atomname (numres , ratm)=atomname 
atomnum(numres , ratm)=i atom 
prvnumres=numres 
end if 
end if 
end do 



return 
end 

function s_eqi (strngl, strng2) 
integer i 
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integer lenl 
integer len2 
integer lenc 
logical s_eqi 
character si 
character s2 

character (len=*) strngl 
character (len=*) strng2 

lenl=len(strngl) 
Ien2=len(strng2) 
1 enc=mi n (1 enl , 1 en2) 

s_eqi=. false. 

do i=l, lenc 

sl=strngl(i :i) 

s2=strng2(i :i) 

call c_cap(sl) 

call c_cap(s2) 

if(sl.ne.s2) then 
return 

end if 
end do 

do i=lenc+l, lenl 

if (strngl(i :i) .ne. " ') then 
return 

end if 
end do 

do i=lenc+l, len2 

if (strng2(i :i) .ne. ' ') then 
return 

end if 
end do 

s_eqi=.true. 

return 
end 



subroutine c_cap(c) 
character c 
integer itemp 

itemp = ichar ( c ) 

i f (97 . 1 e . i temp . and . i temp . 1 e . 122) then 

c=char(itemp-32) 
end if 

return 
end 

c========================================================== 

subrouti ne f i nd_cms (coord , res^atomname , res_atomnum , 
$ numres , maxres , mxratm) 

integer maxres 
integer mxratm 

real*8 coord (max res, mxratm, 3) 

character (1 en=4) res_atomname (maxres , mxratm) 

character (len=4) res_name (maxres) 
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integer res_atomnum(maxres) 



real*8 am(maxres,mxratm) 

integer numres 

real*8 totalm 

real*8 xl 

real*8 x2 

real*8 x3 

do ires=l, numres 

do iatom=l, res_atomnumCi res) 

i f (res_atomname (i res , i atom) (1 : 1) . eq . ' H ' ) 

am(i res , i atom)=1837 . 15dO 
i f ( res_atomname (i res , i atom) (1 : 1) . eq . ' c * ) 

am(i res , i atom) =21874 . 66d0 
i f ( res_atomname (i res , i atom) (1 : 1) . eq . ' N ' ) 

am(ires,iatom)=25526.04dO 
if (res_atomname (ires, i atom) (1:1) .eq. 'O') 

am(i res , i atom)=29156. 94d0 
i f ( res_atomname (i res , i atom) (1 : 1) . eq . ' S ' ) 

am(i res , i atom) =58444 . 168d0 

end do 
end do 

total m=0. do 
xl=0.dO 
x2=0.d0 
x3=0.dO 

do ires=l, numres 

do iatom=l, res_atomnum(i res) 
total m=total m+am(i res , i atom) 
xl=xl+am(i res , i atom) *coord (i res , i atom , 1) 
x2=x2+am (i res , i atom) *coord (i res , i atom , 2) 
x3=x3+am(i res , i atom) *coord (i res , i atom , 3) 
end do 
end do 

xl=xl/totalm 

x2=x2/totalm 
x3=x3/totalm 

pri nt* , ' ***************************** ' 
print*,' Center of Mass ' 

print*, 'x=' ,xl 
print*, 'y=' ,x2 
print*, 'z=' ,x3 

pri nt* , ' ***************************** ■ 

xl=coord(27,17,l) 
x2=coord(27,17,2) 
x3=coord(27,17,3) 

xl=30.dO 
x2=14.d0 
x3=-14.d0 

do ires=l, numres 

do iatom=l, res_atomnum(i res) 

coord (i res , i atom , l)=coord (i res , i atom , 1) -xl 
coord (i res , i atom , 2)=coord (i res , i atom , 2) -x2 

Page 42 



main program.txt 
coord (i res , i atom , 3)=coord (i res , i atom , 3) -x3 
end do 
end do 

return 
end 



subroutine find_charge(res_atomname, res^atomnum, res^name, 
$ charge , endcharge , numres , maxres , mxratm) 

integer maxres 
integer mxratm 

realms coord(maxres, mxratm, 3) 

cha racte r (1 en=4) res_atomname (max res , mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



integer charge (maxres) 
integer endcnarge(2) 
integer numres 
integer i res 
integer iatom 
integer count_H 
integer count^O 
logical chg 

c* initialize the each residue's charge *c 

charge(l:maxres)=0 
chg=.true. 

do ires=l, numres 
count_H=0 
count_O=0 
chg=.true. 

c* for ASP residue to check the charge *c 

if (res_name(i res) .eq. 'ASP') then 
do iatom=l, res_atomnum(i res) 

if(res_atomname(i res, iatom) (1:2) -eq. 'HD') chg=. false, 
end do 

if (chg) charge(i res)=-l 
end if 

c* for GLU residue to check the charge *c 

if(res_name(ire5) .eq. 'GLU') then 
do iatom=l, res_atomnum(i res) f 

if (res_atomname(i res, iatom) (1:2) -eq. 'HE') chg=. false, 
end do 

if (chg) charge(ires)=-l 
end if 

c* for LYS residue to check the charge *c 

if (res_name(i res) .eq. 'LYS') then 
do iatom=l, res_atomnum(i res) 

if (res_atomname (ires, iatom) (1:2) .eq. 'HZ') count_H=count_H+l 
end do 

if (count_H.eq.3) charge(ires)=l 
end if 

c* for ARG residue to check the charge *c 

if (res_name(i res) .eq. 'ARG') then 
do iatom=l, res_atomnum(i res) 

if(res_atomname(ires,iatom)(l:2) .eq. 'HH') count_H=count_H+l 
end do 
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if(count_H.eq.4) charge (i res j=l 
end if 

c* for HIS/HID residue to check the charge *c 

if(res_name(i res) .eq. 'HIS' .or . res_nameO res) .eq. 'HID') then 
do iatom=l, res_atomnum(i res) 

if (res_atomname(i res, i atom) (1:3) .eq. 'HD" .or. 
$ res_atomname(ires,iatom)(l:3) .eq. 'HE') count_H=count_H+l 

end do 

if (count_H.eq.4) charge(i res)=l 
end if 

c* determine the charge for N-end residue in the chain *c 

count_H=0 
count__o=0 

if(ires.eq.l) then 

if (res_name(i res) .ne. 'PRO') then 
do iatom=l, res_atomnum(i res) 

i f (res_atomname(i res , i atom) (1 : 2) 
$ res_atomname (i res , i atom) (1:2) 

$ res_atomname(i res , i atom) (1 : 2) 

$ res_atomname(i res , iatom) (1:2) 

$ res_atomname(i res , i atom) (1 : 2) 

$ res_atomname(i res , i atom) (1 : 2) 

$ res_atomname(i res , i atom) (1 : 2) 

end do 

if (count_H.eq.3) endcharge(l)=l 
else 

do iatom=l, res_atomnum(i res) 

if (res_atomname(i res, iatom) (1: 3) 
$ res_atomname(i res, iatom) (1:3) 

$ res_atomname(i res, iatom) (1:3) .eq 

$ res_atomname(i res , i atom) (1:3). eq 

$ res_atomname (i res , i atom) (1 : 3) . eq 

end do 

if (count„H.eq.2) endcharge(l)=l 
end if 
end if 

c* determine the charge for c-end residue in the chain *c 

if (ires.eq.numres) then 

do iatom=l, res_atomnum(ires) 

if(res_atomname(ires,iatom)(l:3) .eq. 'OXT') then 

endcharge(2)=-l 
end if 
end do 
end if 



.eq 
.eq 
.eq 
.eq 
■eq 
■eq 
.eq 



.eq 
.eq 



'Hi' 
'H2' 
'H3' 
'IH' 
•2H' 
'3H' 



or. 
or. 
or. 
or. 
or. 
or. 



'H') count_H=count_H+l 



•IH' 
'2H' 



or. 
or. 



'Hi" .or. 
'H2' .or. 

'H') count_H=count_H+l 



end do 



return 
end 



subroutine change_coord (coord, res_atomnum,atomnum,numres, 

maxres,mxratm) 



integer maxres 
integer mxratm 



real*8 coord(maxres , mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

character (len=4) res^name (maxres) 

integer res_atomnum(maxres) 
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integer atomnum(maxres,mxratm) 
integer iatom 
integer numres 
integer ires 
character (len=6) wl 
character (len=4) atomname 
integer num(lOOO) 
real x(lOOO) 
real y(lOOO) 
real 2(1000) 
integer i 

open(3333,file='HIvl.pdb' ,status=*old*) 
do iatom=l, 985 

read(3333,"(a6,i5,a3,16x,f8.3,f8.3,f8.3)") 
$ wl, num(i atom) , atomname, X (iatom) , y (i atom) ,z (iatom) 

end do 

do ires=l, numres 

do iatom=l, res_atomnum(i res) 
do i=l, 985 

if (atomnum(i res, iatom) .eq.num(i)) then 
coord (i res , i atom , 1) =x (i ) 
coord (i res, i atom, 2)=y(i) 
coord (i res , i atom, 3)=z(i ) 
end if 
end do 
end do 
end do 

return 
end 

0=====================================================:==========: 

subroutine print_pdb(coord, res_atomname, res_atomnum, 
$ res_name , numres , maxres , mxratm) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

characte r (1 en=4) res^atomname (maxres , mxratm) 

character(len=4) res_name (maxres) 

integer res_atomnum(maxres) 



integer ires 
integer numres 
integer iatom 
real occ 
real temp 
integer atomnum 

occ=1.00 

temp=0.00 

atomnum=0 

do ires=l, numres 

do iatom=l, res_atomnum(i res) 
atomnum=atomnum+l 

i f(res_atomname(i res, i atom) (1 : 1) .eq. '1' .or. 
$ res_atomname(i res, iatom) (1:1) .eq. '2' .or. 

$ res_atomname(ires,iatom)(l:l) .eq. '3') then 
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Write(5000,"(a4,2x,i5,lx,a4,lx,a3,2x,i4,4x,3f8.3,2f6.2)") 
$ • ATOM ' , atomnum , res.atomname (i res , i atom) , res_name(i res) , 

$ i res , coord (i res , i atom ,1:3), occ , temp 

else 

write(5000,"(a4,2x,i5,2x,a4,a3,2x,i4,4x,3f8.3.2f6-2)") 
$ 'ATOM' ,atomnum,res_atomname(i res, i atom) , res_name(i res) , 

$ ires, coord (i res , i atom ,1:3), occ , temp 

end if 
end do 
end do 

write(5000,"(a3,3x,i5,2x,a4,a3,2x,i4)") 'TER' ,atomnum+l, ' ' , 
$ res_name(numres) , numres 

write(5000,"(a3)") 'END' 



c= 
c= 



return 
end 



subrouti ne cal sel ect (Theory , Basi sSet , 1 evel ) 



integer level 

character (len=20) Theory 

character (len=20) Basisset 



wri te(* 
writeC* 
writeC* 
writeC* 
writeC* 
writeC* 
writeC* 
readC*, 
writeC* 
writeC* 
readC*, 
writeC* 
writeC* 



'Protein Cut levels:" 



1. CH3NH2 

2. CH2RNH2 

3. CH2RNHC0H 

4. NH2COCHRNH2 

5. NH2C0CHRNHC0H' 
*) 'Your Choice: ' 
) level 

*) 'YOU have chosen #', level , 'cut method' 
*) 'Select Theory: ' 



0 Theory 

,*) 'You have chosen 
,*) 'select Basisset: 
readC*,*) Basisset 
writeC*,*) 'You have chosen 



Theory 



Basi sSet 



return 
end 



subroutine current_residueCcoord, res_atomname, res_atomnum, 
$ Frag , Frag^atomname , Frag_atomnum , 

$ i res , maxres , mxratm) 

integer maxres 
integer mxratm 

real*8 coord Cmaxres, mxratm, 3) 

character CI en=4) res_atomname Cmaxres, mxratm) 

character CI en=4) res_nameCmaxres) 

integer res_atomnumCmaxres) 



real*8 Frag Cmaxres, mxratm, 3) 

character CTen=4) Frag_atomnameCmaxres, mxratm) 

integer Frag_atomnumCmaxres) 

integer Frag_id Cmaxres) 



integer i 
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integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character(len=4) Al 
character (len=4) a2 
character (len=4) aa 



do iatom=l, res_atomnum(i res) 
do k=l, 3 

F rag (i res , i atom , k) =coord (i res , i atom , k) 
end do 

Frag_atomname(i res , i atom)= 
$ res_atomname (i res , i atom) (1 : 1) 

end do 

Frag_atomnum(i res)=res_atomnum(i res) 

return 
end 

c============================================================== 

subrouti ne CHBnext (coord , res_atomname , res_atomnum , i res , 
$ Frag , Frag_atomname , Frag^atomnum , j res , 

$ maxres , mx ratm , step) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 
character(len=4) res_atomname (max res, mxratm) 
characterOen=4) res_name(maxres) 
integer res_atomnum(maxres) 



real*8 Fraq(maxres, mxratm, 3) 
character(len=4) Frag_atomname(maxres, mxratm) 
integer Frag_atomnum(maxres) 
integer Frag_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character Oen=4) Al 
character (len=4) A2 
character(len=4) AA 



A1='CA' 

A2='C' 

AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
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$ i r es+1 , i res+1 , Al , A2 , max res , mx ratm , x , y) 

call boncnength(x,y,AA) 

A1='CA' 
A2='CB' 

cal 1 twopoi nts (coord , res^atomnum , res_atomname , 
$ i res+1, i res+l,Al, A2,maxres,mxratm,x,yl) 

call bondlength(x,yl,AA) 
step=0 

do iatom=l, res_atomnum(i res+1) 

i f(res_atomname(i res+1, i atom) .eq. 'N' .or. 



$ res_atomname(i res+1, i atom) .eq. 'H' .or. 

$ res_atomname(i res+1, i atom) .eq. 'C .or. 

$ res^atomname (i res+1 , i atom) . eq . ' CA ' . or . 

$ res_atomname (i res+1 , i atom) . eq . ' HA ' . or . 

$ re5_atomname(i res+1, i atom) .eq. 'CB' .or . 

$ res_atomname(i res+1, i atom) .eq. 'ha2' .or. 

$ res_atomname(i res+1, i at om) .eq. 'HA3') then 

step=step+l 
do k=l, 3 

FragC j res , Frag_atomnum( j res)+step , k)= 
$ coord(i res+1, i atom, k) 



end do 

if Cres_atomname(i res+1, iatom) .eq. 'C') then 

Frag_atomname( j res , Frag_atomnum( j res)+step)= ' H ' 

do k=l, 3 

FragCjres, Frag_atomnum(j res)+step, k)=yCk) 
end do 
end if 

if (res_atomnameCi res+1, iatom) .eq- 'CB') then 

Frag_atomname(j res , Frag_atomnum(j res)+step)=' H ' 
do k=l, 3 

FragC j res , Frag_atomnum( j res)+step , k)=yl(k) 
end do 
end if 

if (res_atomname(i res+1, iatom) .ne. 'C' .and. 
$ res_atomname(i res+1, iatom) .ne. 'CB') then 

Frag_atomname(j res , Frag_atomnum(j res)+step)= 
$ res_atomname(i res+1 , i atom) (1:1) 

end if 
end if 
end do 

Frag_atomnum(j res)=Frag_atomnum(j res)+step 

return 
end 

0=========================:===========:=========================== 

subrouti ne CH2Rnext (coord , res_atomname , res^atomnum, i res , 
$ Frag , Frag_atomname , Frag^atomnum, j res , 

$ maxres,mxratm,step) 

integer maxres 
integer mxratm 

real*8 coord (maxres, mxratm, 3) 

character(l en=4) res^atomname (maxres , mxratm) 

characterO en=4) res_name (maxres) 

integer res_atomnum(maxres) 



real*8 Frag(maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnum(maxres) 
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integer Frag_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) aa 



Al='CA' 

A2='C' 

AA='CC' 

call twopoints (coord , res_atomnum, res_atomname, 
$ i res+l,i res+l,Al,A2,maxres,mxratm,x,y) 

cal 1 bondl ength (x , y , AA) 
step=0 

do iatom=l, res_atomnum(i res+1) 

if (res_atomname(i res+1, i atom) . ne. *0' .and . 
$ res_atomname(i res+1, i atom) .ne. 'OXT') then 

step=step+l 
do k=l, 3 

FragCj res , Frag_atomnum( j res)+step, k)= 
$ coo rd (i res+1 , i atom , k) 

end do 

if (res_atomnameCi res+1, i atom) .eq. 'C') then 

Frag_atomname(j res , Frag_atomnum( j res)+step)= ' H ' 
do k=l, 3 

FragCj res , Frag_atomnum(j res)+step , k)=y (k) 
end do 
else 

Frag_atomname(j res , Frag_atomnum( j res)+step)= 
$ res_atomname (i res+1 , i atom) (1:1) 

end if 
end if 
end do 

Frag_atomnum(j res)=Frag_atomnum(j res)+step 

return 
end 

subrouti ne CONHnext (coord , res_atomname , res^atomnum , i res , 
$ Frag , Frag_atomname , Frag_atomnum , j res , 

$ maxres,mxratm,step) 

integer maxres 
integer mxratm 

real*8 coord (max res, mxratm, 3) 

character (1 en=4) res_atomname(maxres , mxratm) 

character (1 en=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Frag (maxres, mxratm, 3) 
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character (len=4) Frag_atomname(maxres,mxratm) 
integer Frag_atomnum(maxres) 
integer Frag_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character Oen=4) AA 



Al='N' 
A2='CA' 
AA= • CN • 

cal 1 twopoi nts (coord , res_atomnum, res_atomname , 
$ i res+2,i res+2,Al,A2,maxres,mxratm,x,y) 

cal 1 bondl ength (x , y , AA) 
step=0 

do iatoni=l, res_atomnum(i res+1) 
step=step+l 
do k=l, 3 

FragCj res , Frag_atomnum(j res)+step, k)= 
$ - coord (i res+1, 1 atom, k) 

end do 

Frag_atomname( j res , Frag_atomnum(j res)+step)= 
$ res_atomname(i res+1, i atom) (1: 1) 

end do 

do iatom=l, res_atomnum(i res+2) 

if(res_atomnameCires+2,iatom).eq. 'N' .or. 
$ res_atomname(i res+2,iatom) -eq. 'H* .or. 

$ res^atomname (i res+2 , i atom) . eq . ' CA ' ) then 

step=step+l 
do k=l, 3 

FragCj res , Frag_atomnum(j res)+step , k)= 
$ coord (i res+2, i atom, k) 

end do 

if (res_atomname(i res+2 ,iatom) .eq. 'CA') then 

Frag_atomname( j res , Frag_atomnum( j res)+step)= ' H ' 
do R=l, 3 

FragCj res , Frag^atomnumCj res)+step , k)=y Ck) 
end do 
else 

Frag^atomnameCj res , Frag_atomnumCj res)+step)= 
$ res_atomnameCi res+2, i atom) CI :1) 

end if 
end if 
end do 

Frag_atomnumC j res)=Frag_atomnumC j res)+step 

return 
end 



subroutine NH2prevCcoord, res_atomname, res_atomnum,i res, 
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$ Frag , Frag^atomname , Frag_atomnum , j res , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (len=4) res_atomname (maxres, mxratm) 

character (len=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Frag (maxres , mxratm, 3) 
character(Ten=4) Frag_atomname(maxres, mxratm) 
integer Frag_atomnum(maxres) 
integer Frag_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
realms y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) a2 
character (len=4) aa 



Al='N" 

A2='c' 

AA='CN' 

cal 1 twopoi nts (coord , res_atomnum, res^atomname, 
$ ires,ires-l,Al,A2, max r es , mx r atm , x , y) 

call bondlength(x,y,AA) 
step=0 

do iatom=l, res_atomnum(i res-1) 

if(res_atomname(ires-l,iatom) .eq, 'C') then 

step=step+l 
do k=l, 3 

Frag(j res , Frag_atomnum(j res)+step , k)=y (k) 
end do 

Frag_atomname( j res , Frag_atomnum( j res)+step)= ' H ' 
end if 
end do 

Frag_atomnum( j res)=Frag_atomnum( j res)+step 

return 
end 

C===================================================:============ 

subroutine NHCOprev(coord, res_atomname, res_atomnum,i res, 
$ Frag , Frag_atomname , Frag_atomnum , j res , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 
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character(l en=4) res_atomname(maxres , mxratm) 
character(len=4) res_name(maxres) 
integer res_atomnum(maxres) 



real*8 Frag(maxres, mxratm, 3) 

character (Ten=4) Frag_atomname(maxres, mxratm) 

integer Frag_atomnum(maxres) 

integer Frag_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) a2 
character (len=4) AA 



A1='C' 

A2='CA' 

AA='CC' 

cal 1 twopoi nts (coord , res_atomnum , res_atomname , 
$ i res-1 , i res-1 , Al , a2 , maxres , mxratm , x , y) 

cal 1 bondl engt h (x , y , AA) 
step=0 

do iatom=l, res_atomnum(i res-1) 

if (res_atomname(i res-1, iatom) .eq. 'C' .or. 
$ res_atomname(i res-1, iatom) .eq. '0' .or. 

$ res_atomname(i res-1, iatom) .eq. 'CA') then 

step=step+l 
do k=l. 3 

FragC j res , Frag_atomnum( j res)+step , k)= 
$ coord (i res-1, iatom, k) 

end do 

if (res_atomname(i res-1, iatom) .eq. 'CA') then 

Frag_atomname( j res , Frag_atomnum(j res)+step) = ' H ' 
do k=l, 3 

FragC j res , Frag_atomnum(j res)+step , k)=y (k) 
end do 
else 

Frag_atomname(j res , Frag_atomnum( j res)+step)= 
$ res_atomname (i res-1 , i atom) (1:1) 

end if 
end if 
end do 

Frag_atomnum(j res)=Frag_atomnum(j res)+step 

return 
end 

c============================================================== 

subroutine Cap_CH3 (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ maxres , mxratm , step) 



integer maxres 
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integer mxratm 

real*8 coord (max res, mxratm, 3) 

character (len=4) res^atomname (max res, mxratm) 

character (1 en=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

character (1 en=4) Cap_atomname (max res , mxratm) 

integer cap_atomnum(maxres) 

integer Cap_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
rea1*8 yl(3) 
character (len=4) Al 
character (len=4) a2 
character (len=4) AA 



Al='CA' 

A2='C' 

AA='CC' 

cal 1 twopoi nts (coo rd , res_atomnum , res^atomname , 
$ i res , i res , Al , a2 , maxres , mxratm , x , y) 

cal 1 bondl ength (x , y , AA) 
Al=*CA' 
A2='CB' 

cal 1 twopoi nts (coord , res_atomnum , res^atomname , 
$ i res , i res , Al , A2 , maxres , mxratm , x , yl) 

call bondl ength(x,yl,AA) 

do iatom=l, res_atomnum(i res) 

if (res_atomname(i res, i atom) .eq. 'N' .or . 
$ res_atomname(i res,iatom) .eq. 'H' .or. 

$ res_atomname(i res, i atom) .eq. 'C* .or . 

$ res_atomname(i res, i atom) .eq. 'CA' .or . 

$ res_atomname(i res, i atom) .eq. 'ha' .or . 

$ res_atomname(i res, i atom). eq. 'CB' .or. 

$ res_atomname(i res,iatom) .eq. 'HA2' .or. 

$ res_atomname (i res , i atom) . eq - ' ha3 ' ) then 

step=step+l 

do k=l, 3 

Cap(i res , step , k)=coord(i res , i atom , k) 
end do 

if (res_atomname(i res , i atom) . eq . 'C') then 
Cap_atomname (i res , step)= ' H ' 
do k=l, 3 

cap(i res , step , k)=y (k) 
end do 

end if 

if (res_atomname(i res, i atom) .eq. 'CB') then 
Cap_atomname (i res , step)= ' H ' 
do k=l, 3 

Cap(i res , step , k)=yl(k) 
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end do 
end if 

if (res_atomname(i res, i atom) .ne. 'c' .and. 
res_atomname (i res , i atom) . ne . ' CB ' ) then 
Cap_atomname(i res , step)= 

res_atomname(i res , iatom) (1:1) 

end if 
end if 
end do 

return 
end 



subroutine Cap_CH2R(coord, res^atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap^atomnum , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 
character(len=4) res^atomname (maxres, mxratm) 
character (1 en=4) res_name(maxres) 
integer res_atomnum(maxres) 



$ 
$ 



real*8 Cap(maxres, mxratm, 3) 

character (len=4) Cap_atomname (maxres, mxratm) 

integer Cap_atomnum(maxres) 

integer Cap_id (maxres) 

integer i 
integer j 
integer iatom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
rea1*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) AA 



Al='CA' 

A2='C' 

AA='CC' 

call twopoints(coord, res_atomnum, res_atomname, 
$ i res,i res, Al,A2, maxres, mxratm, x,y) 

cal 1 bondl ength(x , y , aa) 

do iatom=l, res_atomnum(i res) 

if (res_atomname(i res, iatom) .ne. '0' .and. 
$ res_atomname(i res, iatom) .ne. 'OXT') then 

step=step+l 
do k=l, 3 

cap (i res , step , k)=coord (i res , i atom , k) 
end do 

if (res_atomname(i res, iatom) .eq. 'c') then 
Cap_atomname (i res , step)= ' H ' 
do k=l, 3 
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Cap(i res , step , k)=y (k) 
end do 
else 

Cap_atomname(i res,step)= 
$ res_atomname (i res , i atom) (1:1) 

end if 
end if 
end do 

return 
end 



subroutine Cap_NH2 (coord, res_atomname, res_atomnum, 
$ ires, Cap , Cap_atomname , Cap.atomnum , 

$ maxres , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (1 en=4) res_atomname (maxres , mxratm) 

character (1 en=4) res_name (maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

character (len=4) Cap_atomname (max res, mxratm) 

integer Cap_atomnum(maxres) 

integer Cap_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character Oen=4) Al 
character (len=4) A2 
character (len=4) AA 



A1='N' 
A2='C' 
AA='CN' 

call twopoi nts (coord, res_atomnum, res_atomname, 
$ i res , i res-1 , Al , a2 , maxres , mxratm , x , y) 

call bondlength(x,y,AA) 

do iatom=l, res_atomnum(i res-l) 

i f (res_atomname (i res-1 , i atom) . eq . ' c ' ) then 
step=step+l 
do k=l, 3 

Cap(i res , step , k)=y (k) 
end do 

Cap_atomname (i res , step)= ' H ' 
end if 
end do 

return 
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end 

C====================================:====================== 

subrouti ne Cap_NHCO (coord , res_atomname , res_atomnum, 
$ ires, Cap , Cap_atomname , Cap_atomnum , 

$ max res .mxratm, step) 

integer maxres 
integer mxratm 

real*8 coord(maxres, mxratm, 3) 

character (len=4) res_atomname(maxres, mxratm) 

characte r (1 en=4) res_name (maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

character (len=4) Cap_atomname (max res, mxratm) 

integer Cap_atomnum(maxres) 

integer cap_id(maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
characterOen=4) Al 
character (len=4) A2 
character (len=4) aa 



Al='c' 

A2='CA' 
AA=*CC* 

call twopoi nts (coord, res_atomnum, res_atomname, 
$ i r e s - 1 , i r e s - 1 , Al , a2 , max r e s , mx r atm , x , y ) 

call bondlength(x,y,AA) 

do iatom=l, res_atomnum(i res-1) 

if (res_atomname(i res-1, i atom) .eq. 'c' .or. 
$ res_atomname(i res-1, i atom) .eq. '0' .or. 

$ res_atomname(i res-1, iatom) -eq. 'CA") then 

step=step+l 
do k=l, 3 

Cap (i res , step , k)=coord (i res-1 , i atom , k) 
end do 

if (res_atomname(i res-1, iatom) .eq. 'CA') then 

Cap_atomname (i res , step)= ' H ' 

do k=l, 3 

cap(i res, step, k)=y(k) 

end do 
else 

Cap_atomname(i res , step)= 
$ res_atomname(i res-1, iatom) (1:1) 

end if 
end if 
end do 

return 
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end 



subroutine Cap^CONH (coord, res_atomname, res_atomnum, 
$ i res , Cap , Cap_atomname , Cap_atomnum , 

$ max res , mxratm , step) 

integer maxres 
integer mxratm 

real*8 coordCmaxres, mxratm, 3) 

character (1 en=4) res_atomname(maxres , mxratm) 

character (len=4) res_name (maxres) 

integer res_atomnum(maxres) 



real*8 Cap(maxres, mxratm, 3) 

character (len=4) Cap_atomname (max res, mxratm) 

integer cap_atomnum(maxres) 

integer cap_id (maxres) 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
real*8 yl(3) 
character (len=4) Al 
character (len=4) A2 
character (len=4) AA 



Al='N* 

A2='CA' 

AA='CN' 

call twopoi nts (coord, res_atomnum, res_atomname, 
$ i res+1 , i res+1 , Al , A2 , maxres , mxratm , x , y) 

call bondlength(x,y ,AA) 

do iatom=l, res_atomnum(i res) 
step=step+l 
do k=l, 3 

cap (i res , step , k)=coord (i res , i atom , k) 
end do 

Cap_atomname(i res , step)= 
$ res_atomname (i res , i atom) (1:1) 

end do 

do iatom=l, res_atomnum(i res+1) 

if (res_atomname(i res+1, i atom) .eq. 'N' .or . 
$ res__atomname(i res+1, i atom). eq. 'H' .or. 

$ res_atomname (i res+1, i atom). eq. 'CA') then 

step=step+l 
do k=l, 3 

Cap(i res , step , k)=coord (i res+1 , i atom , k) 
end do 

if (res_atomname(i res+1, i atom) .eq. 'CA') then 
Capiat omname (i res , step) = ' H ' 
do k=l, 3 

Cap(i res , step , k)=y (k) 
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end do 
else 

Cap_atomname(i res,step)= 

res^atomname (i res+1 , i atom) (1 : 1) 

end if 
end if 
end do 

return 
end 



subrouti ne get_pccharge(charge , endcharge , pcharge , ccharge , 

res_name , numres , max res , 1 evel , cut) 

character (len=40) Theory 
nteger level 
nteger error 
nteger cut 
integer numres 
integer ncpu 

nteger maxres 

nteger charqe(maxres) 

nteger endcnarge(2) 

nteger pcharge(maxres) 

nteger ccharge(maxres) 

integer ires 
character (1 en=4) res_name(maxres) 

i f (1 evel . eq . 1) then 
do ires=l, numres 

i f (i res . eq . 1) pcharge(i res)=charge(i res) 
$ +endcharge(l) 

if (i res. eq. numres) pcharge(i res)=charge(i res) 
$ +endcharge(2) 

i f (1 . 1 1 . i res . and . i res .It. numres) pcharge (i res)= 
$ charge(ires) 
end do 
end if 

if (level .eq. 2. or. level .eq. 3. or. level .eq.4.or. 
$ level .eq.5) then 

do ires=l, numres 
cssssssssssssssssss cut=0 i.e. cut CA-N bond ssssssssssssssssssc 
if(cut.eq.O) then 

i f (i res - eq . 1) pcharge(i res)=charge(i res)+ 
$ chargeCi res+l)+enacharge(l) 

i f (i res . eq . numres) pcharge(i res)=charge(i res)+ 
$ endcnarqe(2) 

i f (1 . 1 1 . i res . and . i res . It . numres) pcharge (i res)= 
$ charge(i res)+charge(i res+1) 

end if 

cssssssssssssssssss cut=l i.e. cut CA-C bond ssssssssssssssssssc 
if(cut.eq.l) then 

i f (i res . eq . 1) pcharge(i res)=charge(i res) 
$ +endcharge(l) 

if (ires. eq. numres) pcharge(i res)=charge(i res) 
$ +charge(i res-l)+endcharge(2) 

i f (1 . 1 1 - i res . and . i res .It. numres) pcharge (i res)= 
$ charge(i res)+charge(i res-1) 

end if 

cssssssssssssssssss cut=2 i.e. cut C-N bond sssssssssssssssssssc 
if(cut.eq.2) then 
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if (level .eq. 2) then 

i f (i res . eq . 1) pcharge (i res)=charqe(i res) + 
$ chargeCi res+l)+endcnarge(l) 

if (i res -eq.numres; pcharge(i res)=charge(i res)+ 
$ endcharge(2) 

i f (1 . 1 1 . i res . and . i res .It. numres) pcharge(i res)= 
$ chargeCi res)+charge(i res+1) 

end if 

if (level .eq. 3) then 

i f (i res . eq . 1) pcharge(i res)=charge(i res) + 
$ endcharge(l) 

i f (i res . eq . numres) pcharge(i res)=charge(i res) + 
$ charge(i res-l)+endcharqeC2) 

i f (1 . 1 1 . i res . and . i res .It. numres) pcharge (i res)= 
$ chargeCi res)+charge(i res-l) 

end if 

if (level .eq. 4) then 

if Ci res.eq. 1) pchargeCi res)=charqeCi res)+ 
$ chargeti res+l)+endcnargeCl) 

if Ci res. eq.numres) pchargeCi res)=chargeCi res)+ 
$ chargeCi res-l)+endchargeC2) 

if Cl-lt.i res. and. 1 res. It. numres) pchargeCi res)= 
$ chargeCi res)+chargeCi res+l)+chargeCi res-l) 

end if 
end if 
end do 
end if 

csssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 
do ires=2, numres 

i f CI eve! . eg . 1) then 
cchargeCi res)=0 
else 

if Ccut .eq.O) cchargeCi res)=chargeCi res) 
if Ccut .eq.l) cchargeCi res)=chargeCi res-l) 
ifCcut.eq.2) then 

if Clevel .eq.2) cchargeCi res)=chargeCi res) 
if Clevel .eq.3) cchargeCi res)=chargeCi res-l) 
if (level .eq.4) cchargeCi res)=chargeCi res)+ 
$ chargeCi res-l) 

end if 
end if 
end do 

if Cres_nameCl) .eq. 'PRO') then 
c pchargeC2)=chargeC2)+chargeCl)+endchargeCl) 
c cchargeC2)=cchargeC2)+endchargeCl) 

end if 

if Cres^nameCnumres) .eq. 'pro') then 

pcharge Cnumres-l)=chargeCnumres-l)+chargeCnumres) 
$ +endchargeCnumres) 

cchargeCnumres)=cchargeCnumres)+endchargeCnumres) 
end if 

return 
end 

C================:=============================================== 

subrouti ne combi neuni ts Ccoord , res_atomname , res^atomnum , 
$ Corr,corr_atomname,Corr_atomnum,i res, level , 

$ charge , endcharge , cor recharge , maxres , mxratm , numres) 



integer maxres 
integer mxratm 
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real*8 coord(maxres,mxratm,3) 

character(l en=4) res_atomname(maxres ,mxratm) 

characterO en=4) res_name(maxres) 

integer res_atomnum(maxres) 



real*8 Corr(maxres,mxratm,3) 

cha ract e r (1 en=4) Co r r_atomname (max r es , mx ratm) 

integer Corr_atomnum(maxres) 

integer Corr_charge(maxres) 

character (len=40) Theory 

integer level 

integer error 

integer cut 

integer numres 

integer ncpu 

integer i 

integer q 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

real*8 y(3) 

real*8 yl(3) 

characterO en=4) Al 

characterO en=4) A2 

characterO en=4) AA 



integer charqe(maxres) 
integer endcnarge(2) 

do iatom=l, res_atomnum(i res) 
do k=l, 3 

Cor r (i res , i atom , k) =coord (i res , i atom , k) 
end do 

Corr_atomname(i res , iatom)= 
$ res_atomname(i res, i atom) (1:1) 

end do 

do iatom=l, res_atomnum(i res+1) 
do k=l, 3 

Corr (i res , res_atomnum(i res)+i atom, k)= 
$ coord (i res+l,iatom,k) 

end do 

Corr_atomname(i res , res_atomnum(i res)+i atom)= 
$ res_atomname(i res+1, iatom) (1:1) 

end do 

Corr_atomnum(i res)=res_atomnum(i res)+ 
$ res_atomnum(i res+1) 

if (ires.eq.l) then 
i f (1 evel . eq . 1) then 

cal 1 CH3next (coord , res_atomname , res_atomnum , 
$ i res+1 , Cor r , Cor r_atomname , Corr_atomnum , 

$ i res , max res , mxratm , step) 

end if 

if (level .eq. 2. or. level .eq.3) then 

call CH2Rnext (coord, res_atomname, res_atomnum, 
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$ i res+1 , corn , Corr_atomname , Cor r_atomnum , 

$ i res , max res , mxratm , step) 

end if 

if (level .eq. 4. or .level .eq. 5) then 

call CONHnext (coord, res_atomname, res_atomnum, 
$ i res+1 , Cor r , Corr_atomname , Corr_atomnum , 

$ i res , maxres , mxratm , step) 

end if 
end if 

if (l.lt.i res.and.i res.lt. (numres-1)) then 
i f (1 evel . eq . 1) then 

call NH2prev (coord, res_atomname, res^atomnum, 
$ i res , Cor r , Cor r_atomname , Cor r_atomnum , 

$ i res, maxres, mxratm, step) 

call CH3next (coord, res_atomname, res_atomnum, 
$ i res+1 , Cor r , Corr_atomname , Cor r_atomnum , 

$ i res , maxres , mxratm , step) 

end if 

if (level .eq.2) then 

call NH2prev(coord, res_atomname,res_atomnum, 
$ i res , Corr , Cor r_atomname , Corr_atomnum , 

$ i res , maxres , mxratm , step) 

cal 1 CHZRnext (coord , res_atomname , res_atomnum , 
$ i res+1 , Corr , Cor r_atomname , Cor r_atomnum , 

$ i res , max res , mxratm , step) 

end if 

if (level .eq.3) then 

call NHCOprev (coord, res_atomname, res_atomnum, 
$ i res , Corr , Cor r_atomname , Cor r_atomnum , 

$ i res , maxres , mxratm , step) 

call CH2Rnext (coord, res_atomname, res_atomnum, 
$ i res+1 , Cor r , cor r_atomname , Corr_atomnum , 

$ i res, maxres, mxratm, step) 

end if 

if (level .eq-4) then 

call NH2prev(coord, res_atomname, res_atomnum, 
$ i res , cor r , cor r_atomname , Corr^atomnum , 

$ i res , maxres , mxratm , step) 

if(ires.eq. (numres-2)) then 

call CH2Rnext(coord, res^atomname, res^atomnum, 
$ i res+1 , Cor r , Cor r_atomname , Cor r_atomnum , 

$ i res , maxres , mxratm , step) 

else 

call CONHnext (coord, res_atomname, res_atomnum, 
$ i res+1 , Corr , Cor r_atomname , Cor r_atomnum , 

$ i res , maxres , mxratm , step) 

end if 
end if 

if (level .eq. 5) then 

call NHCOprev (coord, res_atomname, res_atomnum, 
$ i res , Corr , Cor r_atomname , Cor r_atomnum , 

$ i res, maxres, mxratm, step) 

if (i res. eq. (numres-2)) then 

cal 1 CH2Rnext (coord , res_atomname , res_atomnum , 
$ i res+1 , Cor r , Cor r_atomname , Cor r_atomnum , 

$ i res , maxres , mxratm , step) 

else 

call CONHnext (coord, res_atomname, res^atomnum, 
$ i res+1 , Co r r , Co r r_atomname , Cor r_atomnum , 

$ i res, maxres, mxratm, step) 



end if 
end if 
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end if 

if (i res.eq. (numres-1)) then 

i f (1 evel . eq . 1 . or . 1 evel . eq . 2 . or . 1 evel . eq . 4) then 
call NH2prev(coord, res_atomname, res_atomnum, 
$ i res , Corr , Cor r_atomname , Corr_atomnum, 

$ i res , maxres , mxratm , step) 

end if 

i f (1 evel . eq . 3 . or . 1 evel . eq . 5) then 

call NHCOprev(coord, res_atomname, res^atomnum, 
$ i res , Cor r , Cor r_atomname , Corr_atomnum , 

$ i res, maxres, mxratm, step) 

end if 
end if 

if(ires.eq.l) then 
if (level .eq.l) then 

corr_charge(i res)=charge(i res)+charge(i res+1) 
$ +endcharge(l) 
else 

Corr_charge(i res)=charge(i res)+charge(i res+1) 
$ +charge(i res+2)+endchargeCl) 

end if 
end if 

if (i res.eq. (numres-1)) then 

Corr_charge(i res)=charge(i res)+charge(i res+1) 
$ +endcharge(2) 
end if 

if(l.lt.i res.and.ires.lt. (numres-1)) then 
i f (1 evel . eq . 1) then 

corr_charge(i res)=charge(i res)+charge(i res+1) 
else 

Corr_charge(i res)=charge(i res)+charge(i res+1) 
$ +charge(ires+2) 
end if 
end if 

return 
end 



subroutine distance(x,y,dis) 

real*8 x(3) 
real*8 y(3) 
real*8 dis 
integer k 

do k=l, 3 

dis=dis+(x(k)-y(k))**2 
end do 

dis=dsqrt(dis) 

return 
end 



subrouti ne mi ndi s (tmp , tmp_atomnum , tmp_atomname , i res , numres , 
$ 1 i gand , 1 i gand_atomnum , maxsi ze , maxres , 

$ mxratm,min) 



integer maxres 
integer mxratm 
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real*8 tmp(maxres,mxratm, 3) 
character(len=4) tmp_atomname(maxrGS,mxratm) 
integer tmp_atomnum(maxres) 
integer tmp_id(maxres) 
integer maxsize 

character(len=4) atom_symbol (maxsize) 
real*8 1 i gand (maxsize, 3 j 
integer 1 i gand_atomnum 
integer ligand^charge 

integer i 
integer j 
integer latom 
integer k 
integer step 
integer ires 
integer jres 
real*8 x(3) 
real*8 y(3) 
rea1*8 yl(3) 
character (len=4) Al 
character (len=4) a2 
character (len=4) AA 



integer numres 
real*8 dis 
rea1*8 min 

min=10000.dO 

do iatom=l, tmp_atomnum(i res) 
do j=l, 1igand_atomnum 
do k=l, 3 

x(k)=ligand(j ,k) 
end do 
do k=l, 3 

y(k)=tmp(i res , iatom, k) 
end do 

call distance(x,y,dis) 
if (min, ge. dis) then 

mi n=ais 
end if 
end do 
end do 

return 
end 

0============================================================== 

subrouti ne sel ectg roup (A_i d , i res , nopol ar , pol ar , charged , 
$ min,maxres,id) 

integer maxres 

integer ires 

integer A_id(maxres) 

integer id 

real*8 min 

real nopol ar 

real polar 

real charged 
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id=l 

if (A_id(i res) .eq.l) then 

if (min.gt.nopolar) id=0 
end if 

if (A^id(i res) .eq. -1) then 
if (min.gt.polar) id=0 
end if 

if (A_id(ires) -eq.O) then 

if (mi n.gt. charged) id=0 
end if 

return 
end 

c========================================================== 

subrouti ne pol ari ty (res_name , Frag_i d , Cap_i d , numres , 
$ maxres , 1 evel , cut) 

integer maxres 
character (len=40) Theory 
integer level 
integer error 
integer cut 
integer numres 
integer ncpu 

character(len=4) res_name(maxres) 
integer Frag^id (maxres) 
integer Cap_id(maxres) 
integer id(maxres) 
integer ires 

do ires=l, numres 



if (res_name(i res) .eq. 'ALA' .or. 

$ res_name(i res) .eq. 'VAL' .or. 

$ res_name(ires) -eq. 'leu' .or. 

$ res_name(i res) .eq. 'ILE' .or. 

$ res_name(ires) .eq. 'PRO' .or. 

$ res_name(i res) .eq. 'MET' .or. 

$ res„name(i res) .eq. 'PHE' .or. 

$ res_name(ires) .eq. 'TRP') id(ires)=l 

if(res_name(ires) .eq. 'GLY' .or. 

$ res_name(ires) .eq. 'SER' .or. 

$ res_name(i res) .eq. 'THR' .or. 

$ res_name(ires) .eq. 'CYS' .or. 

$ res_name(ires) .eq. 'CYX' .or. 

$ res_name(ires) .eq. 'ASN' .or. 

$ res_name(i res) . eq . ' GLN ' . or . 

$ res_name(ires) .eq. 'TYR') id(ires)=-l 

if (res_name(i res) .eq. 'ASP' .or. 

$ res_name(i res) . eq . ' GLU ' . or . 

$ res_name(i res) .eq. ' LYS' .or. 

$ res_name(i res) .eq. 'arg' .or. 

$ res_name(i res) .eq. 'HIS' .or. 

$ res_name(ires) .eq. 'HID') id(ires)=0 



end do 



id(l)=0 
id(numres)=0 
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cssssssssssssssssss cut=0 i.e. cut CA-N bond ssssssssssssssssssc 
if(cut.eq.O) then 

if (level .eq.l) then 
do ires=l, numres 

Frag_icl(i res)=id(i res) 
Cap_id(i res)=l 
end do 
else 

do ires=l, numres-1 

i f (i d(i res) . eq . -1 . and - i d(i res+1) . eq . -1) then 

Frag_id(ires)=-1 
else 

Frag_id(i res)=i d(i res)*id(i res+1) 
end if 

cap_id(i res+l)=id(i res+1) 
end do 

Frag_id(numres)=id(numres) 
end if 
end if 

cssssssssssssssssss cut=l i.e. cut CA-C bond ssssssssssssssssssc 
if (cut. eq.l) then 

i f (1 evel . eq . 1) then 
do ires=l, numres 

Frag_i d(i res)=id(i res) 
Cap_id(i res)=l 
end do 
else 

do ires=2, numres 

i f (i d (i res) . eq . -1 . and . i d(i res-1) . eq - -1) then 

Frag_id(i res)=-l 
else 

Frag_i d (i res) =i d (i res) *i d (i res- 1) 
end if 

Cap_id(i res)=id(i res-1) 
end do 

Frag_id(l)=id(l) 
end if 
end if 

cssssssssssssssssss cut=2 i.e. cut C-N bond sssssssssssssssssssc 
if(cut.eq.2) then 

i f (1 evel . eq . 1) then 
do ires=l, numres 

Frag_id(i res)=id(i res) 
Cap_id(i res)=l 
end do 
end if 

if (level .eq.2) then 
do ires=l, numres-1 

if (id(ires) .eq--l.and.id(ires+l) .eq.-l) then 

Frag_id(i res)=-l 
else 

Frag_id(i res)=id(i res)*id(i res+1) 
end if 

Cap_id(i res+l)=id(i res+1) 
end do 

Frag_id(numres)=id (numres) 
end if 

if (level .eq.3) then 
do ires=2, numres 

i f (i d (i res) , eq . -1 . and - i d(i res-l) . eq . -1) then 

Frag_id(i res)=-l 
else 

Frag_id(i res)=id(i res)*id(i res-1) 
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end if 

Cap_id(i res)=id(i res-1) 
end do 

Frag«id(l)=id(l) 
end if 

if (level .eq. 4) then 
do ires=2, numres-1 

i f (i d(i res) . eq . -1 . and . i d(i res-l) - eq . -1) then 

Frag_id(ires)=-1 
else 

Frag_i d (i res)=i d (i res-1) *i d (i res) *i d (i res+1) 
end if 

cap«id(i res)=id(i res-l)*id(i res) 
end do 

Frag_id(l)=id(l)*id(2) 
Frag_id(numres)=id(numres)*id (numres-1) 
end if 

end if 

return 
end 



subroutine bondlenQth(xl,x2,id) 
implicit real*8(a-n,o-z) 
dimension xl(3) ,x2(3),y(2,3),z(2,3) 
character(len=4) id 

Pl=dacos(-l.dO) 

bond=0.dO 

do i=l, 3 

bond=bond+ (xl(i ) -x2 (i ) ) **2 
end do 

bond=dsqrt(bond) 

do j=l, 3 

y(2,j)=x2(j)-xl(j) 

y(l,j)=0.dO 
end do 

theta=y(2,3)/bond 

if (theta.ge.l.dO) then 

theta=l.dO 

end if 

t heta=dacos (t heta) 

phi=abs(y(2,l))/bond/dsin(theta) 
if (phi .ge.l.dO) then 
phi =1. do 
end if 

phi=dacos(phi) 

if(y(2,l) .qt.0.d0.and,y(2,2) .gt.O.dO) phi=phi 
if(y(2,l) . it.0.d0.and.y(2,2) .gt.O.dO) phi=PI-phi 
i f (y (2 , 1) . 1 1 . 0 . do . and . y (2 , 2) . 1 1 . 0 . dO) phi =Pl+phi 
i f (y (2 , 1) . gt . 0 . do . and . y (2 , 2) , 1 1 . 0 . dO) phi=2 . dO*PI-phi 

if(id.eq. 'CO bond=1.09d0 
if(id.eq. 'CN') bond=1.01d0 

y (2 , l)=bond*dsi n(theta) *dcos (phi ) 
y (2 , 2)=bond*dsi n (theta) *dsi n (phi ) 
y(2 , 3)=bond*dcos(theta) 
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do j=l, 3 

x2(j)=y(2.j)+xl(j) 
end do 

return 
end 



subrouti ne twopoi nts (tmp , tmp_atomnum , tmp^atomname , i res , 
$ j res , atomi , atomj , maxres , mxratm , x , y) 

integer maxres 
integer mxratm 

real*8 tmp(maxres , mxratm, 3) 

character (1 en=4) tmp_atomname (maxres , mxratm) 

integer tmp_atomnum(maxres) 

integer tmp_id (max res) 

integer i 

integer j 

integer latom 

integer k 

integer step 

integer ires 

integer jres 

real*8 x(3) 

realms y(3) 

real*8 yl(3) 

character (len=4) Al 

character (len=4) A2 

character (len=4) AA 



character (len=4) atomi 
character (len=4) atomj 

do iatom=l, tmp_atomnum(i res) 

i f (tmp_atomname (i res , i atom) . eq . atomi ) then 
do k=l. 3 

X (k) =tmp (i res , i atom , k) 
end do 
end if 
end do 

do iatom=l, tmp_atomnum(jres) 

i f (tmp_atomname( j res , i atom) . eq . atomj) then 
do k=l, 3 

y (k)=tmp( j res , i atom , k) 
end do 
end if 
end do 

return 
end 
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